"""
Импортируем библиотеки для расчетов.
numpy - для работы с массивами и подготовки данных
matplotlib - для построения графиков
scipy - для решения линейных уравнений
"""
import numpy as np
import matplotlib.pyplot as plt
import scipyЛекция 1. Введение и простые решения
версия 0.6 от 16.09.2025
Хабибуллин Ринат
1 Введение в дисциплину. Моделирование при исследовании скважин и пластов
1.1 Тезисы для лекции
- Управление месторождениями строится на принятии решений. Принятие решений связано с необратимой тратой ресурсов с целью оптимизации некоторой целевой функции. Чаще всего целевая функция связана с экономикой проекта.
- Исследования - инструмент для получения информации для принятия решений. Полезно разделить понятия данные и информация. Данные - то, что записано, информации - данные обработанные в контексте принятия решений.
- Исследования основаны на моделях объекта управления. Модель - иструмент получения информации из данных. Также модель - инструмент принятия решений на основе информации. Модель - упрощенное представление реальности.
- Модели используемые для принятия решений должны адекватно отражать реальность. Обеспечивается настройкой - адаптацией модели на известные данные и информацию. Такой процесс также можно считать исследованием. Почему?
- Принятие решений на основе информации часто осуществляется на основе прогнозов - альтернативных вариантах принятия решений. Модель - инструмент для построения прогнозов.
- В узком смысле - исследование скважины, в частности гидродинамические исследования - соответствуют адаптации специализированной модели - решения уравнения фильтрации. Адаптация готовит модель для построения прогнозов и принятия решений.
- Понимание всей цепочки создания ценности от данных до принятия решений оказывается полезным, например для обоснования необходимости проведения исследований и оценки ценности информации получаемой в ходе исследований
- Кроме понимания всей цепочки важно развитие навыков проведения расчетов - адаптации различных моделей, построения прогнозов с использованием этих моделей. Этому будет посвящена значительная часть курса.
1.2 Инструменты для моделирования
Работа с моделями в ходе курса может проводиться с использованием различного ПО
python для построения аналитических моделей
Для работы можно использовать jupyter notebook или google colab или любой другой сервис для работы со скважинами. Потребуется освоить базовые библиотеки python - numpy, scipy, matplotlib, mpmath, pandas. Также будут использованы некоторые другие специализированные библиотеки.
Макросы MS Excel
Часть расчетов можно выполнить с использованием макросов MS Excel, например Unifloc VBA (https://github.com/unifloc/unifloc_vba)
Гидродинамический симулятор для построения численных моделей
Базовый вариант - тНавигатор. Версия начиная от 19. Можно будет использовать другие симуляторы - РН-КИМ, Schlumberger eclipse, OPM и другие.
Специализированное ПО для анализа ГДИС
Возможно, что то появится в ходе курса.
2 Уравнение фильтрации
Уравнение фильтрации многофазного потока в пористой среде - основная математическая модель пласта используемая в инженерных приложениях. Используется не само уравнение, а разнообразные частные решения, описывающие те или иные ситуации. Наиболее широко известным примером решения уравнения фильтрации, возможно является, численное решение которое строится трехмерными гидродинамическими симуляторами (eclipse, tnavigator и так далее).
Как правило, при выводе уравнения фильтрации используются следующие соотношения:
- закон Дарси - определяет линейную зависимость скорости фильтрации от давления
\[ u_r=-\frac{k}{\mu}\frac{dp}{dr} \tag{2.1}\]
уравнение неразрывности - или закон сохранения массы, показывает что в потоке не могут образовываться разрывы сплошности: \[ \frac{1}{r}\frac{\partial\left(r\rho u_r\right)}{\partial r}=-\varphi\frac{\partial p}{\partial t} \tag{2.2}\]
уравнение состояния - задает зависимость плотности флюида от давления: \[ c_0=\frac{1}{\rho}\frac{\partial\rho}{\partial p} \tag{2.3}\]
Опираясь на эти соотношения уравнение фильтрации в радиальной форме можно привести к виду для величин СИ.
\[ \frac{\partial ^2 p }{\partial r^2} + \frac{1}{r} \frac{\partial p}{\partial r} = \frac{\varphi \mu c_t}{k} \frac{\partial p}{\partial t} \tag{2.4}\]
Здесь используются следующие обозначения:
- \(u_r\) - скорость фильтрации в направлении \(r\), м/сек
- \(k\) - проницаемость, м\(^2\)
- \(\mu\) - вязкость флюида, Па с
- \(p\) - давление, Па
- \(r\) - расстояние, м
- \(\rho\) - плотность флюида, кг/м\(^3\)
- \(\varphi\) - пористость породы, доли единиц.
- \(c_t\) - общая сжимаемость породы и флюида, 1/Па
- \(t\) - время, сек
Иногда бывает удобно записать уравнение используя другие единицы измерения, например практические метрические - уравнение примет вид:
Уравнение фильтрации в практических метрических единицах измерений
\[ \frac{\partial ^2 p }{\partial r^2} + \frac{1}{r} \frac{\partial p}{\partial r} = \frac{\varphi \mu c_t}{0.00036 k} \frac{\partial p}{\partial t} \tag{2.5}\]
- \(k\) - проницаемость, мД
- \(\mu\) - вязкость флюида, сП
- \(p\) - давление, атм
- \(r\) - расстояние, м
- \(\rho\) - плотность флюида, кг/м\(^3\)
- \(\varphi\) - пористость породы, доли единиц.
- \(c_t\) - общая сжимаемость породы и флюида, 1/атм
- \(t\) - время, час
Вывод уравнения фильтрации основан на следующих предположениях:
- пласт однородный и изотропный - пористость и проницаемость одинаковы во всем пласте и во всех направлениях и не зависят от давления
- добывающая скважина вскрывает весь продуктивный горизонт и обеспечивает радиальный приток к скважине
- пласт насыщен одним флюидом на всем протяжении
- температура не меняется в пласте, изотермичность
2.1 Вывод уравнений фильтрации
2.1.1 Вывод уравнения фильтрации для линейного потока
Уравнение фильтрации основано на трех основных уравнениях:
- законе сохранения массы
- зависимости связывающей градиент давления с расходом, в простейшем виде это закон Дарси
- уравнении состояния - зависимости связывающей сжимаемость флюида (и отчасти среды) и давлением
Выпишем эти уравнения, чтобы разобраться. В уравнениях будем придерживаться обозначений принятых в международной литературе, в частности в публикациях общества инженеров нефтяников SPE. Описание обозначений, а также единиц измерения можно найти в Эрлагер (2006). Особое внимание на данном этапе стоит обратить на размерные коэффициенты в уравнениях. С одной стороны записывая уравнения в СИ их можно упростить и избежать части сложностей с переводными коэффициентами. С другой стороны во многих источниках переводные коэффициенты интесивно используются и это оказывается удобно при программной реализации алгоритмов. В данном описании приведем перевод коэффициентов в практические метрические единицы.
Эрлагер, Роберт. 2006. Гидродинамические методы исследования скважин. М.–Ижевск: Институт компьютерных исследований.
2.1.1.1 Закон сохранения массы
В дифференциальной форме, для потока в направлении \(x\)
\[ -f \dfrac{\partial \rho q_x}{\partial x} = A \dfrac{\partial(\rho \phi)}{\partial t} \tag{2.6}\]
где
- \(f\) - размерный коэффициент, для практических метрических единиц \(f= 0.04167\), для американских промысловых \(f = 0.23394\)
- \(\rho\) - плотность флюида, кг/м\(^3\), lbm/ft\(^3\)
- \(q_x\) - объемный расход в направлении \(x\), м\(^3\)/сут, bbl/day, учитываем для практических метрических единиц \(1\) м\(^3\)/сут \(= 0.04167\) м\(^3\)/час и для промысловых американских \(1\) bbl/day \(= 0.23394\) ft\(^3\)/hr
- x - координата, м, ft
- \(\phi\) - пористость, доли единиц
- \(t\) - время, часы
- \(A\) - площадь, м\(^2\), ft\(^2\)
В перечне указаны размерности для промысловых метрических и американских промысловых единиц измерений, там где они отличаются.
Для практических метрических единиц (2.6) можно записать в виде \[ -0.04167 \dfrac{ \partial (\rho q_x) }{\partial x} = A \dfrac{\partial(\rho \phi)}{\partial t} \tag{2.7}\]
2.1.1.2 Закон Дарси
Простейшим вариантом задания зависимости между измением давления и потоком флюида является линейная зависимость - закон Дарси. При этом линейный множитель в зависимости - функция как коллектора (проницаемость), так и свойств флюида (вязкость).
Закон Дарси в направлении \(x\). \[ q_x = - \dfrac{k_x A}{f\mu} \dfrac{\partial p}{\partial x} \tag{2.8}\]
где
- \(f\) - размерный коэффициент, для практических метрических единиц \(f= 115.74\), для американских промысловых \(f = 887.2\)
- \(k_x\) - проницаемость в направлении движения потока, мД
- \(A\) - площадь, м\(^2\), ft\(^2\)
- \(\mu\) - динамическая вязкость, сП
- \(p\) - давление, атм, psi
- \(x\) - координата движения потока, м
В практических метрических единицах измерения (2.8) будет иметь вид \[ q_x = - \dfrac{k_x A}{115.74\mu} \dfrac{\partial p}{\partial x} \tag{2.9}\]
Пытаясь проверить корректность вывода переводных коэффициентов необходимо учесть что для единиц измерение СИ получим \(f=1\), а необходимые переводные коэффициенты можно записать как:
- \(q\): \(1\) [м\(^3\)/сек] = \(543439\) [bbl/day] = \(86400\) [м\(^3\)/сут]
- \(k\): \(1\) [м\(^2\)] = \(1.01325 \cdot 10^{15}\) [мД]
- \(\mu\): \(1\) [Па\(\cdot\)с] = \(1000\) [сП]
- \(p\): \(1\) [Па] = \(0.0001450\) [psi] = \(0.00000987\) [атм]
- \(x\): \(1\) [м] = \(3.28\) [ft]
- \(A\): \(1\) [м\(^2\)] = \(10.76\) [ft\(^2\)]
2.1.1.3 Уравнение состояния
В простейшем варианте, который мы будет изучать, предполагается, что сжимаемость постоянна и мала
\[ c_t = c_{formation} + c_{fluid} = const \tag{2.10}\]
- \(c_t\) - общая сжимаемость (total compressibility)
- \(c_{formation}\) - сжимаемость породы (formation compressibility), \(c_{formation} = \dfrac{1}{\phi}\dfrac{\partial \phi}{\partial p}\)
- \(c_{fluid}\) - сжимаемость флюида (fluid compressibility), \(c_{fluid} = \dfrac{1}{\rho}\dfrac{\partial \rho}{\partial p}\)
2.1.1.4 Уравнение фильтрации
Подставим выражение для закона Дарси (2.9) в уравнение мат баланса (2.7)
\[ 0.04167 \dfrac{ \partial}{\partial x} { \left[ \rho \dfrac{k_x A}{115.74\mu} \dfrac{\partial p}{\partial x} \right] } = A \dfrac{\partial(\rho \phi)}{\partial t} \tag{2.11}\]
Упрощая (2.5) получим
\[ \frac{\partial (\rho \phi)}{\partial t} = 0.00036 k_x \dfrac{ \partial}{\partial x}\left[ \dfrac{\rho}{\mu}\dfrac{\partial p}{\partial x} \right] \tag{2.12}\]
Уравнение (2.12) подходит для потока газа. Для жидкости, предполагая малую сжимаемость флюида можно его упростить до вида
\[ \frac{\partial p}{\partial t} = 0.00036 \dfrac{k_x}{\phi \mu c_t} \dfrac{ \partial^2 p}{\partial x^2} \tag{2.13}\]
Напомним, здесь
- \(p\) - давление, атм
- \(t\) - время, час
- \(k_x\) - проницаемость в направлении движения потока, мД
- \(\mu\) - динамическая вязкость, сП
- \(\phi\) - пористость, д.е.
- \(c_t\) - сжимаемость, 1/атм
- \(x\) - расстояние, м
Приведенное уравнение (2.13) соответствует линейному потоку в пласте. Практический интерес представляет решение для радиального потока
2.1.2 Вывод уравнения фильтрации для радиального потока
2.1.2.1 Закон сохранения массы
Закон сохранения массы или принцип неразрывности можно выразить в радиальной форме следующим соотношением
\[ \frac{1}{r}\frac{\partial\left(r\rho u_r\right)}{\partial r}=-\varphi\frac{\partial p}{\partial t} \tag{2.14}\]
Принцип неразрывности показывает, что для определенного объема пласта, масса флюида которое втекла в контрольный объем пласта минус масса которая вытекла равна массе которая накопилась в объеме.
Сохранение импульса или закон Дарси можно выразить соотношением
\[ u_r=-\frac{k}{\mu}\frac{d^p}{dr} \tag{2.15}\]
Закон Дарси здесь используется как псевдоустановившаяся аппроксимация обобщенного уравнения сохранения импульса (то есть слагаемым отвечающим за накопление импульса пренебрегаем). Это предположение справедливо если не учитывать возмущения давления в среде двигающиеся со скоростью звука. Все изменения давления описываемые моделью связаны с локальными изменениями градиента давления за счет ламинарного режима потока. Хотя далее в модели будет учтена сжимаемость системы всеми “звуковыми” эффектами в системе мы пренебрегаем.
Поток предполагается горизонтальным, поэтому давление \(p\) может быть использовано в качестве потенциала потока, гравитационными силами пренебрегаем.
Комбинируя уравнения (2.14) , (2.15) получим:
\[ \frac{1}{r}\frac{\partial\left( \dfrac{r\rho k}{\mu}\dfrac{\partial p}{dr}\right)}{\partial r}=\varphi\frac{\partial \rho}{\partial t} \tag{2.16}\]
Уравнение (2.16) – дифференциальное уравнение в частных производных описывающее нестационарный поток однофазного флюида в пористой среде при ламинарном потоке.
Вообще говоря приведенное уравнение является нелинейным, так как плотность \(\rho = \rho(p)\) и вязкость \(\mu = \mu(p)\) являются функциями давления. Уравнение содержит две зависимые переменные - давление \(p\) и плотность \(\rho\). Поэтому для его решения необходимо задать еще одно соотношение, каковым может быть уравнение состояния флюида связывающее плотность флюида и давление \(\rho = \rho(p)\).
2.1.2.2 Флюид постоянной сжимаемости
Нестационарное поведение давления в пласте связано со сжимаемостью системы. При изменении давления в какой то точке, часть флюида сжимается, происходит накопление или отдача флюида, что вызывает задержку в распространении изменения флюида. Несмотря на то, что сжимаемость флюидов и породы малы и во многих случаях ими можно пренебречь, это не верно для пластовых систем для добычи нефти. Большие объемы пласта и флюидов и высокие давления компенсируют малость сжимаемости и требуют ее учета.
Для однофазного флюида разумным является предположение постоянства сжимаемости.
Сжимаемость можно определить как
\[ c=-\frac{1}{V} \left( \frac{ \partial V}{ \partial p} \right) \]
учтем, что
\[ \rho = \frac{m}{V} \]
тогда получим
\[ c=\frac{1}{\rho} \left( \frac{ \partial \rho}{ \partial p} \right) \]
Для флюида с постоянной сжимаемостью, проинтегрировав приведенное уравнение можно получить
\[ \rho = \rho_i e^{c(p-p_i)} \]
где \(\rho_i\) плотность флюида при некотором заданном давлении \(p_i\)
Продифференцировав выражение для плотности по времени получим
\[ c \rho \frac{\partial p}{\partial t} = \frac{\partial \rho}{\partial t} \]
Подставив это выражение в ранее полученное уравнение фильтрации (2.16) получим
\[ \frac{1}{r}\frac{\partial\left( \dfrac{r\rho k}{\mu}\dfrac{\partial p}{dr}\right)}{\partial r}=\varphi c \rho \frac{\partial p}{\partial t} \tag{2.17}\]
Приведенное дифференциальное уравнение в частных производных все еще нелинейно, поскольку зависит от плотности \(\rho\)
2.1.2.3 Общая сжимаемость
Если пористость не является постоянной величиной и меняется с давлением, тогда слагаемое отвечающее за накопление флюида в пласте можно выразить как
\[ \frac{\partial \varphi \rho}{\partial t} = \varphi \frac{\partial \rho}{\partial t}+ \rho \frac{\partial \varphi }{\partial t} = \varphi c_l \rho \frac{\partial p}{\partial t} + \rho \frac{\partial \varphi }{\partial t} \]
где \(c_l\) сжимаемость жидкости.
Определим сжимаемость породы как
\[ c_f = \frac{1}{\varphi} \frac{\partial \varphi}{\partial p} \]
тогда
\[ \frac{\partial \varphi \rho}{\partial t} = \varphi \rho (c_l + c_f) \frac{\partial p}{\partial t} \]
хотя пористость здесь является функцией давления - в первом приближении мы можем считать ее константой равной пористости при некотором среднем давлении в пласте. Это справедливо для маленькой сжимаемости породы, что верно почти всегда.
Уравнение для сжимаемости можно еще уточнить, учтя что в пласте могут находится различные флюиды - вода и нефть с насыщенностями \(s_w\) и \(s_o\)
тогда
\[ c_l = s_o c_o + s_w c_w \]
тогда можно ввести общую сжимаемость системы
\[ c_t = c_l + c_f = s_o c_o + s_w c_w + c_f \]
Заметим, что проницаемость k в законе Дарси это не абсолютная проницаемость, но относительная проницаемость по нефти при насыщенности водой соответствующей связанной воде.
\[ k= k_o (s_{wc}) \]
2.1.2.4 Линеаризация уравнения фильтрации
Раскрыв производную в левой части (2.17) и предположив, что \(\dfrac{\partial p}{\partial r}\) мало а следовательно слагаемым \(r \rho c_t \left( \dfrac{\partial p}{\partial r} \right)^2\) можно пренебречь, получим линеаризованное уравнения фильтрации
\[ \frac{\partial^2 p}{\partial r^2} + \frac{1}{r} \frac{\partial p}{\partial r}= \frac{\varphi \mu c_t}{k} \frac{\partial p}{\partial t} \tag{2.18}\]
2.2 Безразмерные переменные
Часто для анализа уравнений неустановившейся фильтрации используются безразмерные переменные. Мы будем использовать переменные в виде:
\[ r_D = \frac{r}{r_w} \]
\[ t_D = \frac{kt}{\varphi \mu c_t r_w^2} \]
\[ p_D = \frac{2 \pi kh}{q_s B \mu} \left( p_i - p_{wf} \right) \]
\[ q_D = \frac{q}{q_{ref}} \]
Здесь использованы единицы измерения СИ.
- \(q_s\) - дебит скважины на поверхности, приведенный к нормальным условиям м3/с
- \(\varphi\) - пористость, доли единиц
- \(\mu\) - вязкость нефти в пласте, Па с
- \(B\) - объемный коэффициент нефти, м3/м3
- \(p_i\) - начальное давление в пласте, Па
- \(p_{wf}\) - давление забойное, Па
- \(c_t\) - общая сжимаемость системы в пласте, 1/Па
Использование безразмерных переменных позволяет упростить уравнение фильтрации, которое примет вид
\[ \frac{\partial p_D}{ \partial t_D} = \frac{1}{r_D} \frac{ \partial{ \left( r_D \dfrac{\partial p_D}{ \partial r_D} \right) } }{ \partial{r_D} } \]
\[ \frac{\partial p_D}{ \partial t_D} = \dfrac{1}{r_D} \left[ \dfrac{\partial}{\partial r_D} \left( r_D \dfrac{ \partial p_D} {\partial r_D} \right) \right] \]
Решение этого уравнения - функция безразмерного давления от безразмерных времени и расстояния \(p_D(r_D, t_D)\)
Для практических расчетов удобнее бывает использовать безразмерные переменные полученные для практических метрических единиц измерения.
Определение безразмерных переменных в практических метрических единицах
\[ r_D = \frac{r}{r_w} \]
\[ t_D = \frac{0.00036 kt}{\varphi \mu c_t r_w^2} \]
\[ p_D = \frac{kh}{ 18.42 q_s B \mu} \left( p_i - p_{wf} \right) \]
\[ q_D = \frac{q}{q_{ref}} \]
Здесь использованы практические метрические единицы измерения.
- \(q_s\) - дебит скважины на поверхности, приведенный к нормальным условиям м3/сут
- \(\varphi\) - пористость, доли единиц
- \(\mu\) - вязкость нефти в пласте, сП
- \(B\) - объемный коэффициент нефти, м3/м3
- \(p_i\) - начальное давление в пласте, атм
- \(p_{wf}\) - давление забойное, атм
- \(c_t\) - общая сжимаемость системы в пласте, 1/атм
Уравнение фильтрации для радиального потока в линеаризованном виде можно записать в виде
\[ \frac{\partial p}{\partial t} = 0.00036 \dfrac{k}{\phi \mu c_t} \dfrac{1}{r} \left[ \dfrac{\partial}{\partial r} \left( r \dfrac{ \partial p} {\partial r} \right) \right] \]
Использование безразмерных переменных позволяет упростить уравнение фильтрации, которое примет вид
\[ \frac{\partial p_D}{ \partial t_D} = \dfrac{1}{r_D} \left[ \dfrac{\partial}{\partial r_D} \left( r_D \dfrac{ \partial p_D} {\partial r_D} \right) \right] \]
Решение этого уравнения - функция безразмерного давления от безразмерных времени и расстояния \(p_D(r_D, t_D)\)
Стационарное решение в безразмерных переменных примет вид \[ p_D = \ln r_{eD} - \ln r_D \]
2.2.1 Расчет безразмерных переменных в Unifloc VBA
Несмотря на простоту определений безразмерных переменных их часто приходится применять при проведении расчетов. Поэтому в надстройке Unifloc VBA реализован набор функций расчета безразмерных переменных.
Эти функции начинаются с префикса transient_def.
transient_def_cd
transient_def_cs_1atm
transient_def_td
transient_def_t_day
transient_def_pd
transient_def_pwf_atmaОписания функций и из аргументов можно найти в руководстве пользователя Unifloc VBA
3 Стационарные решения уравнения фильтрации
Широкое распространение на практике получили стационарные решения уравнения фильтрации. Приведем некоторые из них.
3.1 Решение для постоянного давления на круговой границе
Рассматривается самая простая модель работы добывающей скважины - радиальная стационарная фильтрация в однородном изотропном пласте круговой формы. Скважина находится в центре пласта (Рисунок 3.1). На границе пласта поддерживается постоянное давление. Фактически это означает, что через границу пласта идет поток жидкости, уравновешивающий дебит скважины.
Решение можно получить как решения уравнения фильтрации, учитывая стационарность потока
\[ \frac{\partial ^2 p }{\partial r^2} + \frac{1}{r} \frac{\partial p}{\partial r} = 0 \tag{3.1}\]
Или можно получить его непосредственно из закона Дарси, который должен быть приведен к радиальной форме и в таком варианте известен как Формула Дюпюи. Для приведенной конфигурации можно записать закон Дарси в форме.
\[ u_r=\frac{q}{2\pi rh}=\frac{k}{\mu}\frac{dP}{dr} \]
Проинтегрировав выражение по замкнутому контуру радиуса \(r_e\) вокруг скважины получим выражение известное как формула Дюпюи
\[ q=\frac{2\pi kh\left(P_e-P_w\right)}{\mu\left(\ln{\dfrac{r_e}{r_w}}\right)} \tag{3.2}\]
В приведенном выражения использованы единицы СИ.
Здесь
- \(u_r\) - приведенная скорость фильтрации на расстоянии \(r\) от скважины, м/с
- \(q\) - объемные дебит скважины в рабочих условиях, м\(^3\)/с
- \(r\) - радиус - расстояние от центра скважины, м
- \(r_e\) - радиус зоны дренирования, на котором поддерживается постоянное давление, м
- \(r_w\) - радиус скважины, на котором замеряется забойное давление, м
- \(P\) - давление, Па
- \(P_e\) - давление на внешнем контуре дренирования, Па
- \(P_w\) - давление на забое скважины, Па
- \(k\) - проницаемость, м\(^2\)
- \(\mu\) - вязкость нефти в зоне дренирования, Па с
На практике часто бывает удобнее пользоваться значениями в практических метрических единицах измерения.
\[ q=\frac{kh\left(P_e-P_w\right)}{ 18.4 \mu\left(\ln{\dfrac{r_e}{r_w}}\right)} \tag{3.3}\]
где
- \(q\) - объемные дебит скважины в рабочих условиях, м\(^3\)/сут
- \(r\) - радиус - расстояние от центра скважины, м
- \(r_e\) - радиус зоны дренирования, на котором поддерживается постоянное давление, м
- \(r_w\) - радиус скважины, на котором замеряется забойное давление, м
- \(P\) - давление, атм
- \(P_e\) - давление на внешнем контуре дренирования, атм
- \(P_w\) - давление на забое скважины, атм
- \(k\) - проницаемость, мД
- \(\mu\) - вязкость нефти в зоне дренирования, сП
Далее если не указано особо будем использовать практические метрические единицы.
3.2 Расчет решения с использованием python
Для работы с решениями уравнения фильтрации удобно использовать язык программирования python. Расчет на python может быть реализован на основе библиотек numpy scipy.
Для удобства дальнейшего изложения и использования расчетных функций при создании функций и переменных на языке python названия формируются по следующим принципам:
- сначала указывается, что расчитывается в функции, в данном случае - давление
pилиdp - потом указываются пояснения - в данном случае
dp_ss- steady state pressure - в конце указывается размерность в которой ожидается получаение ответа - в данном случае
atma- абсолютные атмосферы.
"""
Определим функции для расчета стационарного решения
"""
def dp_ss_atm(q_liq_sm3day = 50,
mu_cP = 1,
b_m3m3 = 1.2,
kh_mDm = 40,
r_e_m = 240,
r_m = 0.1):
"""
функция расчета перепада давления в произвольной точке пласта
на расстоянии r_m от центра скважины для стационарного решения
- q_liq_sm3day - дебит жидкости на поверхности в стандартных условиях
- mu_cP - вязкость нефти (в пластовых условиях)
- B_m3m3 - объемный коэффициент нефти
- kh_mDm - kh пласта
- r_e_m - радиус контрура питания, м
- r_m - расстояние на котором проводится расчет, м
"""
return 18.42 * q_liq_sm3day * mu_cP * b_m3m3/ kh_mDm * np.log(r_e_m/r_m)
def p_ss_atma(p_res_atma = 250,
q_liq_sm3day = 50,
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10,
r_e_m = 240,
r_m = 0.1):
"""
функция расчета давления в произвольной точке пласта
на расстоянии r_m от центра скважины для стационарного решения
- p_res_atma - пластовое давление, давление на контуре питания
- q_liq_sm3day - дебит жидкости на поверхности в стандартных условиях
- mu_cP - вязкость нефти (в пластовых условиях)
- B_m3m3 - объемный коэффициент нефти
- k_mD - проницаемость пласта
- h_m - мощность пласта, м
- r_e_m - радиус контрура питания, м
- r_m - расстояние на котором проводится расчет, м
"""
return p_res_atma - dp_ss_atm(q_liq_sm3day = q_liq_sm3day,
mu_cP = mu_cP,
b_m3m3 = b_m3m3,
kh_mDm = k_mD * h_m,
r_e_m = r_e_m,
r_m = r_m)Функции расчетов могут быть использованы для построения графиков, например с использованием matplotlib
Показать код
"""
Построим график распределения давления в пласте
"""
# формируем массив расстояний для которых будем проводить расчет
r_arr = np.linspace(0.1, 100, 500)
# рассчитываем массив давлений на соответствующих расстояниях
# для расчета используется векторный расчет numpy - нет необходимости делать цикл в явном виде
# для примера показана передача всех аргументов созданной функции
p_arr = p_ss_atma(p_res_atma = 250,
q_liq_sm3day = 50,
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10,
r_e_m = 240,
r_m = r_arr)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(5,8))
# рисуем график в обычных координатах
ax1.plot(r_arr, p_arr) # команда отрисовки графика по заданным массивам
ax1.plot(-r_arr, p_arr) # отрицательная ветка
# настраиваем график
ax1.set_xlabel('r, m')
ax1.set_ylabel('p, atma')
# рисуем график в логарифмических координатах
ax2.plot(r_arr, p_arr) # команда отрисовки графика по заданным массивам
ax2.plot(-r_arr, p_arr) # отрицательная ветка
# настраиваем график
ax2.set_xlabel('r, m')
ax2.set_xscale('symlog', linthresh=0.1, linscale=0.6)
plt.show()
3.3 Учет скин-фактора для стационарного решения
Скин-фактор — гидродинамический параметр, характеризующий дополнительное фильтрационное сопротивление течению флюидов в околоскважинной зоне пласта, приводящее к изменению добычи (дебита) по сравнению с совершенной (идеальной) скважиной. Скин-фактор может приводить как к снижению дебита (например при загрязнении ПЗС), так и увеличению (образование высокопроводящих каналов в ПЗС).
Концепция скин-фактора получила широкое распространение на практике. Все инженеры-нефтяники знают этот параметр и оперируют им на практике.
Изначально скин-фактор был введен как параметр учитывающий изменение проницаемости (загрязнение) призабойной зоны при расчете производительности скважины. Такое загрязнение может быть вызвано различными причинами:
- проникновением бурового раствора в пласт и блокировкой поровых каналов;
- набуханием глин при контакте с фильтратом бурового раствора;
- химическим осаждением элементов бурового раствора, жидкости глушения или пластовых флюидов в призабойной зоне скважины, например осаждением солей или асфальтенов;
- продвижением песчаных частиц к стволу скважины;
- повреждением породы при перфорации;
- другими причинами.
Скин-фактор определяется как безразмерный перепад давления на стенке скважины или в призабойной зоне пласта.
\[ S = \Delta p_{skin} \dfrac{kh}{18.41 q_{sf} \mu} \]
- \(S\) - скин
- \(\Delta p_{skin}\) - перепад давления за счет скин-фактора
- \(k\) - проницаемость, мД
- \(h\) - эффективная мощность пласта, м
- \(q_{sf}\) - дебит на забое, м\(^3\)/сут
- \(\mu\) - вязкость, сП
Такая формулировка позволяет записать стационарное уравнение притока к скважине в виде
\[ p_D = \ln{r_e{D}} - \ln{r_D} +S \]
Скин-фактор во всех расчетах можно заменить на эквивалентный радиус скважины. Идея заключается в том, что изменение проницаемости в призабойной зоне скважины можно представить как измененный радиус скважины в пласте постоянной проницаемости. Такой вариант работает и для положительных и для отрицательных значений скин-фактора, что иногда бывает удобно с вычислительной точки зрения \[ S = -\ln \left(\dfrac{r_{eff.w}}{r_w} \right) \] или \[ r_{eff.w} = r_w e^{-S} \]
- \(r_{eff.w}\) - Эффективный радиус скважины,м
- \(r_{w}\) - радиус скважины,м
- \(S\) - скин
Для модели загрязненной призабойной зоны величину скин-фактора можно выразить формулой Хокинса (Hawkins 1956). Скин-фактор для плоскорадиального установившегося потока несжимаемой жидкости:
Hawkins, Murray F. 1956. «A Note on the Skin Effect». Journal of Petroleum Technology 8 (12): 65–66. https://doi.org/10.2118/732-G.
\[ S =\left( \frac{k}{k_s} -1\right)\ ln\frac{r_s}{r_w} \tag{3.4}\]
здесь:
- \(k_s\) - проницаемость в загрязненной ПЗП;
- \(k\) - однородная проницаемость по всему пласту;
- \(r_s\) - радиус загрязненной зоны;
- \(r_w\) - радиус скважины.
Концепция скин-фактора оказалась удобной для описания характеристики соединения скважины и пласта и была распространена на другие случаи, когда производительность скважины могла отличаться от производительности идеальной скважины:
- для горизонтальных скважин;
- для скважин вскрывающих пласт под углом;
- для скважин пересеченных трещиной ГРП;
- для скважин вскрытых перфорацией и учета гидравлического сопротивления потока на перфорационных отверстиях;
- другими причинами.
Для многих подобных случаев предположение о радиальном притоке к скважине не верно, но величину скин-фактора используют, так как она позволяет сравнить производительность скважины со сложным заканчиванием с простой вертикальной скважиной. В таких случая говорят о псевдорадиальном скин-факторе - такой величине скин-фактора \(S\), которая обеспечила бы такую же производительность для вертикальной скважины полностью вскрывающей пласт.
Для стационарной радиальной модели притока учет скин-фактора приведен к следующим соотношениям:
\[ (P_e - P_{wf}) = \frac{18.42\mu q }{\ k h}(\ln\frac{r_e}{r_w}+S) \tag{3.5}\]
\[ q=\frac{kh\left(P_e-P_w\right)}{ 18.42 \mu\left(\ln{\dfrac{r_e}{r_w}} + S\right)} \tag{3.6}\]
3.4 Решение для постоянного давления на круговой границе с учетом среднего давления в области дренирования
Показать код
# оценим среднее давление в области дренирования
r_external_m = 240
r_well_m = 0.1
p_average_atma = scipy.integrate.quad(lambda r:p_ss_atma(p_res_atma=250,
q_liq_sm3day=50,
mu_cP=1,
b_m3m3=1.2,
k_mD=40,
h_m=10,
r_e_m=r_external_m,
r_m=r) * 2 * r,
r_well_m,
r_external_m)[0] / r_external_m**2
print(f"p_average_atma = {p_average_atma}")
r_arr = np.linspace(r_well_m, r_external_m, 500)
r_full = np.linspace(-r_external_m, r_external_m, 10)
# рассчитываем массив давлений на соответствующих расстояниях
# для расчета используется векторный расчет numpy - нет необходимости делать цикл в явном виде
# для примера показана передача всех аргументов созданной функции
p_arr = p_ss_atma(p_res_atma=250,
q_liq_sm3day=50,
mu_cP=1,
b_m3m3=1.2,
k_mD=40,
h_m=10,
r_e_m=r_external_m,
r_m=r_arr)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8,8))
# рисуем график в обычных координатах
ax1.plot(r_arr, p_arr) # команда отрисовки графика по заданным массивам
ax1.plot(-r_arr, p_arr)
ax1.plot(r_full, np.full_like(r_full, p_average_atma), linestyle='dashed')
# настраиваем график
#ax1.set_title('Воронка депрессии - распределение давления в пласте')
#ax1.set_xlabel('r, m')
ax1.set_ylabel('p, atma')
# рисуем график в логарифмических координатах
ax2.plot(r_arr, p_arr) # команда отрисовки графика по заданным массивам
ax2.plot(-r_arr, p_arr)
ax2.plot(r_full, np.full_like(r_full, p_average_atma), linestyle='dashed')
# настраиваем график
#ax2.set_title('Воронка депрессии - распределение давления в пласте')
ax2.set_xlabel('r, m')
ax2.set_ylabel('p, atma')
ax2.set_xscale('symlog', linthresh=0.1, linscale=0.6)
plt.show()p_average_atma = 248.6184605726843
В приведенное выражение входит значение давления на контуре, которым не всегда бывает удобно пользоваться. В практических случаях значение на контуре трудно оценить, контур зоны дренирования может значительно отличаться от кругового, да и радиус оценить может быть сложно. Удобнее пользоваться средним давлением в зоне дренирования \(\bar{P}\), которое может быть оценено по материальному балансу (смотри рисунок 3.3).
Вывод уравнения фильтрации для постоянного давления на границе с использованием среднего давления
В этом случае выражение для дебита примет вид \[ q=\frac{kh\left( \bar{P}-P_w\right)}{ 18.42 \mu\left(\ln{\dfrac{r_e}{r_w}} - \dfrac{1}{2}+ S \right)} \tag{3.7}\]
3.5 Решение для круговой непроницаемой границы
Схема модели радиального притока для условия непротекания на круговой границе приведена на рисунке 3.4
При условии непротекания давления на границе условия стационарности (неизменности давления) не достигаются. При работе скважины с постоянным дебитом забойное давление будет постоянно снижаться. Однако начиная с некоторого момента, когда влияние скважины достигнет границ - давление в всей области дренирования начнет снижаться равномерно (смотри Рисунок 3.5).
Такой режим, при котором забойное давление меняется, но перепад давления \(P_e - P_w\) остается постоянным называют псевдо-установившимся режимом работы (pss - pseudo steady state).
Для псевдо-установившегося режима можно записать выражение
\[ q=\frac{kh\left(P_e-P_w\right)}{ 18.42 \mu\left(\ln{\dfrac{r_e}{r_w}} - \dfrac{1}{2} + S\right)} \tag{3.8}\] где
- \(q\) - объемные дебит скважины в рабочих условиях, м3/сут;
- \(r\) - радиус - расстояние от центра скважины, м;
- \(r_e\) - радиус зоны дренирования, на котором поддерживается постоянное давление, м;
- \(r_w\) - радиус скважины, на котором замеряется забойное давление, м;
- \(P\) - давление, атм;
- \(P_e\) - давление на внешнем контуре дренирования, атм;
- \(P_w\) - давление на забое скважины, атм;
- \(k\) - проницаемость, мД;
- \(\mu\) - вязкость нефти в зоне дренирования, сП.
3.6 Решение для круговой непроницаемой границы с учетом среднего давления в зоне дренирования
Аналогично случаю для постоянного давления на границе можно переписать выражение с использованием среднего давления в области дренирования.
\[ q=\frac{kh\left( \bar{P}-P_w\right)}{ 18.42 \mu\left(\ln{\dfrac{r_e}{r_w}} - \dfrac{3}{4}+ S \right)} \tag{3.9}\]
Вывод уравнений для псевдо-установившегося режима работы
3.7 Стационарные решения для вертикальной скважины в резервуаре произвольной формы
Здесь уравнения и методы расчета для горизонтальных, наклонно направленных скважин, скважин с ГРП, горизонтальных скважин с МГРП.
\[ q=\frac{kh\left( \bar{P}-P_w\right)}{ 18.42 \mu\left(\ln{\dfrac{2.2458 A}{C_A r_w^2}} + S \right)} \tag{3.10}\]
- \(q\) - объемные дебит скважины в рабочих условиях, м3/сут;
- \(A\) - площадь области дренирования, м\(^2\);
- \(C_A\) - фактор формы, зависит от формы резервуара и расположения скважины;
- \(r_w\) - радиус скважины, на котором замеряется забойное давление, м;
- \(\bar{P}\) - среднее давление в области дренирования, атм;
- \(P_e\) - давление на внешнем контуре дренирования, атм;
- \(P_w\) - давление на забое скважины, атм;
- \(k\) - проницаемость, мД;
- \(\mu\) - вязкость нефти в зоне дренирования, сП.
%todo -надо переделать рисунок и добавить ссылки на источники
Безразмерное время достижения псевдо-установившегося режима притока, определяемое видом резервуара
\[ t_{DApss} = \frac{kt_{pss}}{\varphi \mu c_t A} \tag{3.11}\]
3.8 Формула Дюпюи
Простое решение для задачи стационарного притока к вертикальной скважине в однородном изотропном пласте круговой формы с постоянным давлением на границе имеет вид
\[ Q=\dfrac{kh}{18.42\mu B} \dfrac{P_{res}-P_{wf}}{\ln \dfrac{r_e}{r_w} + S} \]
где:
- \(Q\) - дебит скважины на поверхности, приведённый к нормальным условиям, ст. м\(^3\)/сут
- \(\mu\) - вязкость нефти в пласте, сП
- \(B\) - объёмный коэффициент нефти, м\(^3\)/м\(^3\)
- \(P_{res}\) - пластовое давление или давление на контуре с радиусом \(r_e\), атма
- \(P_{wf}\) - давление забойное, атма
- \(k\) - проницаемость, мД
- \(h\) - мощность пласта, м
- \(r_e\) - внешний контур дренирования скважины, м
- \(r_w\) - радиус скважины, м
- \(S\) - скин-фактор скважины, безразм.
Это решение известно как закон Дарси или формула Дюпюи.
Выражение можно переписать в виде
\[ P_{r} = P_{res} - 18.42\dfrac{ Q\mu B }{kh} \left[ \ln\dfrac{r_e}{r} +S \right] \tag{3.12}\]
который удобен для расчёта распределения давления в пласте \(P_r\) на произвольном расстоянии от скважины \(r\). В выражении (2) задано граничное значение давления \(p_e\) на контуре \(r_e\). Расчёт позволит найти любое значение внутри контура, в том числе и забойное давление \(P_{wf}\) на \(r=r_w\)
Выражение можно переписать \[ P_{r} = P_{wf} + 18.42\dfrac{ Q\mu B }{kh} \left[ \ln\dfrac{r}{r_w} +S \right] \tag{3.13}\]
где по известному дебиту и забойному давлению можно найти давление в пласте. При известном пластовом давлении можно оценить радиус контура на котором оно достигается.
3.8.1 Формула Дюпюи в декартовых координатах
Для построения карты распределения давлений в пласте полезно вспомнить, что расстояние от скважины с координатами \((x_{well}, y_{well})\) до произвольной точки пласта с координатами \((x,y)\) можно найти по формуле
Показать код
"""
Построим карту или сетку распределения давления в пласте
с использованием векторных функций numpy
"""
# зададим параметры воронки депрессии
pres = 250
r_e = 300
# зададим координатную сетку основываясь на параметрах
x = np.linspace(-r_e, r_e, 300)
y = np.linspace(-r_e, r_e, 300)
# рассчитаем вспомогательные вектора для построения сетки
xv, yv = np.meshgrid(x, y)
# зададим координаты скважины
xwell1 = 0
ywell1 = 0
# рассчитаем значение давлений во всех точках сетки
# расчет ведется для матрицы координата с использованием векторных возможностей numpy
p_mesh = p_ss_atma(r_m=((xv-xwell1)**2 + (yv-ywell1)**2)**0.5, p_res_atma=pres, r_e_m=r_e)
# удалим значения за контуром, так как в данном случае они не имеют смысла
p_mesh[np.where(p_mesh > pres)] = pres
# построим отображение в виде контурной карты
plt.rcParams["figure.figsize"] = (7,7)
plt.contour(x, y, p_mesh, levels=100)
plt.show()
\[ r=\sqrt{ (x-x_{well})^2 + (y-y_{well})^2 } \tag{3.14}\]
Тогда выражение для расчета давления в любой точке пласта примет вид
\[ P_{r} = P_{res} - 18.42\dfrac{ Q\mu B }{kh} \left[ \ln\dfrac{r_e}{\sqrt{ (x-x_{well})^2 + (y-y_{well})^2 }} +S \right] \tag{3.15}\]
Простой вариант расчета - можно создать пустую матрицу со значениями давления по сетке и перебирая все точки на сетке/матрице рассчитать давления
3.9 Суперпозиция для стационарного решения
3.9.1 Суперпозиция для нескольких скважин с постоянным дебитом
Показать код
"""
Построим карту или сетку распределения давления в пласте
для ряда/галереи скважин
"""
# зададим параметры воронки депрессии
pres = 250
r_e = 300
# зададим координаты скважины (всего будет 5 скважин у которых меняется только х координата)
xwell = [-100, -50, 0 , 50, 100]
ywell = 0
# зададим координатную сетку основываясь на параметрах
x = np.linspace(-r_e, r_e, 300)
y = np.linspace(-r_e, r_e, 300)
# рассчитаем вспомогательные вектора для построения сетки
xv, yv = np.meshgrid(x, y)
# зададим пустой список матриц с перепадами давлений от каждой отдельной скважины
p_mesh_i=[]
# для каждой скважины найдем ее влияние на давления
for xi in xwell:
# рассчитаем матрицу расстояний от элементов сетки до скважины i
r_well = ((xv-xi)**2 + (yv-ywell)**2)**0.5
# рассчитаем значение давлений во всех точках сетки для скважины i
p_mesh_i_ = p_ss_atma(r_m=r_well, p_res_atma=pres, r_e_m=r_e)
# удалим значения за контуром, так как в данном случае они не имеют смысла
p_mesh_i_[np.where(p_mesh_i_ > pres)] = pres
# сохраним влияние скважины i в списке матриц влияния отдельных скважин
p_mesh_i.append(p_mesh_i_)
# найдем суммарное влияние все скважин
p_mesh_sum = 0
for p_mesh_i_ in p_mesh_i:
# найдем сумму решений, помним что суммировать можно депрессии
p_mesh_sum = pres - p_mesh_i_ + p_mesh_sum
p_mesh_sum = pres - p_mesh_sum
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4,8))
# построим отображение в виде контурной карты
ax1.contourf(x, y, p_mesh_sum, levels=100)
# построим отображение в виде контурной карты
ax2.contour(x, y, p_mesh_sum, levels=15)
plt.show()
Для стационарного решения работает принцип суперпозиции - сумма двух решений также будет решением, это позволяет построить карту для нескольких скважин. Давление в любой точке пласта можно найти по формуле
\[ P_{res} - P_{x,y} = \sum_{i} 18.42\dfrac{ Q_i\mu B }{kh} \left[ \ln\frac{r_e}{\sqrt{ (x-x_{w.i})^2 + (y-y_{w.i})^2 }} +S \right] \tag{3.16}\]
Выражение (3.16) справедливо только если \(\sqrt{ (x-x_{w.i})^2 + (y-y_{w.i})^2 }< r_e\).
Расчет распределения давления в пласте показан на рисунке 3.8. Интересно отметить, что при постоянном дебите нескольких скважин на одной линии забойное давление на них будет отличаться, что ножно наглядно показать построив распределение давления по линии расположения скважин (Рисунок 3.9).
Показать код
"""
Построим график изменения давления по линии расположения скважин
"""
xline = np.linspace(-300, 300, 10000)
yline = 0
pline = np.zeros_like(xline)
for xi in xwell:
r_well = ((xline-xi)**2 + (yline-ywell)**2)**0.5
r_well = np.where(r_well < 0.1, 0.1, r_well)
pline = pline + p_ss_atma(r_m=r_well,
p_res_atma=0,
r_e_m=r_e,
q_liq_sm3day = 50,
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10)
pline = pres + pline
plt.plot(xline, pline)
plt.show()
3.9.2 Суперпозиция для нескольких скважин с постоянным забойным давлением
При наличии нескольких скважин можно записать выражение для оценки забойных давлений скважин
Показать код
# зададим параметры воронки депрессии
pres = 250
r_e = 300
# зададим координатную сетку основываясь на параметрах
x = np.linspace(-r_e, r_e, 500)
y = np.linspace(-r_e, r_e, 500)
# рассчитаем вспомогательные вектора для построения сетки
xv, yv = np.meshgrid(x, y)
# зададим координаты скважины (всего будет 5 скважин у которых меняется только х координата)
xwell = [-100, -50, 0 , 50, 100]
ywell = [0,0,0,0,0]
pwell = [100, 100, 100, 100, 100] # забойные давления скважин
qwell = [0,0,0,0,0] # дебиты скважин, должны быть найдены
# создадим заготовки матрицы А и вектора В по формулам (4) и (5)
A = np.zeros((5,5))
B = np.zeros(5)
# сформируем матрицу А
for i,v in enumerate(xwell):
for j,v in enumerate(ywell):
# найдем расстояния от одной скважины до другой
# ищем расстояния между центрами скважин
r_ij = ((xwell[i]-xwell[j])**2 + (ywell[i]-ywell[j])**2)**0.5
# если расстояние 0 значит ищем влияние скважины саму на себя
# тогда подставляем радиус скважины
if r_ij == 0:
r_ij = 0.1
# чтобы воспользоваться ранее заданной формулой в виде функции
# вызовем ее с нулем пластовым давлением и единичным дебитом
A[i,j] = - p_ss_atma(p_res_atma = 0,
q_liq_sm3day = 1,
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10,
r_e_m = 240,
r_m = r_ij)
B[i] = pres - pwell[i]
# найдем решение = значения дебитов при которых забойные равны заданным значениям
qwell = scipy.linalg.solve(A,B)
# напечатаем найденные дебиты
for i, q in enumerate(qwell):
print(f"q_{i} = {q}")
# зададим пустой список матриц с перепадами давлений от каждой отдельной скважины
p_mesh_i=[]
# для каждой скважины найдем ее влияние на давления
for i,xi in enumerate(xwell):
# рассчитаем матрицу расстояний от элементов сетки до скважины i
r_well = ((xv-xwell[i])**2 + (yv-ywell[i])**2)**0.5
# для красоты отрисовки карты давлений пренебрежем значениями в радиусе одного метра от скважин
r_well[r_well<1]=0.1
# рассчитаем значение давлений во всех точках сетки для скважины i
p_mesh_i_ = p_ss_atma(r_m=r_well, p_res_atma=pres, r_e_m=r_e, q_liq_sm3day=qwell[i])
# удалим значения за контуром, так как в данном случае они не имеют смысла
p_mesh_i_[np.where(p_mesh_i_ > pres)] = pres
# сохраним влияние скважины i в списке матриц влияния отдельных скважин
p_mesh_i.append(p_mesh_i_)
# найдем суммарное влияние все скважин
p_mesh_sum_new = 0
for p_mesh_i_ in p_mesh_i:
# найдем сумму решений, помним что суммировать можно депрессии
p_mesh_sum_new = pres - p_mesh_i_ + p_mesh_sum_new
p_mesh_sum_new = pres - p_mesh_sum_new
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(5,9))
# построим отображение в виде контурной карты
ax1.contourf(x, y, p_mesh_sum_new, levels=100)
# построим отображение в виде контурной карты
ax2.contour(x, y, p_mesh_sum_new, levels=15)
plt.show()q_0 = 263.5244312064029
q_1 = 214.68400868202497
q_2 = 202.93807721505908
q_3 = 214.68400868202494
q_4 = 263.5244312064028
\[ P_{res} - P_{wf.j} = \sum_{i} 18.42\dfrac{ Q_i\mu B }{kh} \left[ \ln\dfrac{r_e}{\sqrt{ (x_{w.j}-x_{w.i})^2 + (y_{w.j}-y_{w.i})^2 }} +S \right] \tag{3.17}\]
Если считать забойные давления \(P_{wf.j}\) известными а дебиты скважин \(Q_i\) не известными, тогда выражение (6) можно рассматривать как систему линейных алгебраических уравнений вида
\[ AX = B \tag{3.18}\]
Где \[ A_{[i,j]} = 18.42\dfrac{ \mu B }{kh} \left[ \ln\dfrac{r_e}{\sqrt{ (x_{w.j}-x_{w.i})^2 + (y_{w.j}-y_{w.i})^2 }} +S \right] \tag{3.19}\]
\[ B_{[j]}=P_{res} - P_{wf.j} \tag{3.20}\]
такую систему можно решить например с использованием пакета scipy.linalg.
Убедиться, что забойные давления на скважинах равны можно построив график распределения давления в пласте вдоль линии скважин (Рисунок 3.11). При этом дебиты на скважинах будут разные. Величины можно увидеть на рисунке 3.10.
Показать код
"""
Построим график изменения давления по линии расположения скважин
"""
xline = np.linspace(-300, 300, 10000)
yline = 0
pline = np.zeros_like(xline)
for i,xi in enumerate(xwell):
r_well = ((xline-xi)**2 + (yline-ywell[i])**2)**0.5
r_well = np.where(r_well < 0.1, 0.1, r_well)
pline = pline + p_ss_atma(r_m=r_well,
p_res_atma=0,
r_e_m=r_e,
q_liq_sm3day=qwell[i],
mu_cP = 1,
b_m3m3 = 1.2,
k_mD = 40,
h_m = 10)
pline = pres + pline
plt.plot(xline, pline)
plt.show()
3.9.3 Расчет поля давлений для нескольких скважин на стационарном режиме с использованием классов
Классы в python помогают задавать сложные структуры данных (инкапсуляция) и работать с ними как с едиными объектами. Скважины с их параметрами, а также группы скважин могут быть заданы с использованием классов.
Класс для расчета влияния одной скважины - объединяет все необходимые данные для скважины - позволяет рассчитать депрессию и давления в любой точке пласта с учетом радиуса влияния и радиуса скважины
Показать код
class Well:
"""
класс Well описывает влияние одной скважины на поле давления
используется стационарное решение
"""
def __init__(self, name, xw_m, yw_m) -> None:
self.x_m = xw_m
self.y_m = yw_m
self.name = name
self.kh_mDm = 10.0
self.b_m3m3 = 1.0
self.mu_cP = 1.0
self.skin = 0
self.r_e_m = 250.0
self.r_w_m = 0.1
self.q_liq_sm3day = 10
def calc_dp_ss_atm(self, x_m, y_m):
"""
расчет перепада давления в произвольной точке пласта
с заданными координатами
для стационарного решения уравнения фильтрации
x_m, y_m координаты в которых рассчитывается давление
"""
r_m = np.sqrt( (self.x_m - x_m)**2 + (self.y_m - y_m)**2 )
r_m = np.where(r_m >= self.r_e_m, self.r_e_m ,r_m)
r_m = np.where(r_m < self.r_w_m, self.r_w_m ,r_m)
return dp_ss_atm(q_liq_sm3day=self.q_liq_sm3day,
mu_cP=self.mu_cP,
b_m3m3=self.b_m3m3,
kh_mDm=self.kh_mDm,
r_e_m=self.r_e_m,
r_m=r_m)
def calc_dp_ss_well_atm(self, w):
return self.calc_dp_ss_atm(w.x_m, w.y_m)
def calc_p_ss_atma(self, x_m, y_m, p_res_atma):
"""
функция расчета давления в произвольной точке пласта
с заданными координатами
для стационарного решения уравнения фильтрации
x_m, y_m координаты в которых рассчитывается давление
"""
return p_res_atma - self.calc_dp_ss_atm(x_m, y_m)
class Wells:
"""
класс для группы взамодействующих скважин
"""
def __init__(self, name:list, xw_m:list, yw_m:list,
kh_mDm=10, b_m3m3=1,
mu_cP=1, r_e_m =250, r_w_m=0.1) -> None:
self.well_dict: dict[str, Well] = {}
for (n, x, y) in zip(name, xw_m, yw_m):
self.well_dict.update({n: Well(n,x,y)})
self.set_kh_mDm(kh_mDm)
self.set_pvt(b_m3m3=b_m3m3, mu_cP=mu_cP)
self.set_r_e_m(r_e_m)
self.set_r_w_m(r_w_m)
self.p_res_atma = 250
def set_qliq_sm3day(self, q_liq_list_sm3day:list):
for (wname, q) in zip(self.well_dict, q_liq_list_sm3day):
self.well_dict[wname].q_liq_sm3day = q
def set_kh_mDm(self, kh_mDm):
for wn in self.well_dict:
self.well_dict[wn].kh_mDm = kh_mDm
def set_r_e_m(self, r_e_m):
for wn in self.well_dict:
self.well_dict[wn].r_e_m = r_e_m
def set_r_w_m(self, r_w_m):
for wn in self.well_dict:
self.well_dict[wn].r_w_m = r_w_m
def set_pvt(self, mu_cP, b_m3m3):
for wn in self.well_dict:
self.well_dict[wn].b_m3m3 = b_m3m3
self.well_dict[wn].u_cP = mu_cP
def calc_p_ss_atma(self, x_m, y_m):
dp = 0
for wn in self.well_dict:
dp = dp + self.well_dict[wn].calc_dp_ss_atm(x_m, y_m)
return self.p_res_atma - dp
def estimate_qliq_from_pwf(self, wn_list:list, pwf_list:list):
dim = len(wn_list) # размерность матрицы для расчета дебитов
# создадим заготовки матрицы А и вектора В по формулам (4) и (5)
A = np.zeros((dim,dim))
B = np.zeros(dim)
for i, wni in enumerate(wn_list):
self.well_dict[wni].q_liq_sm3day = 1
for j, wnj in enumerate(wn_list):
A[i,j] = self.well_dict[wni].calc_dp_ss_well_atm(self.well_dict[wnj])
B[i] = self.p_res_atma - pwf_list[i]
q_list = scipy.linalg.solve(A,B)
for i, wni in enumerate(wn_list):
self.well_dict[wni].q_liq_sm3day = q_list[i]
xlist = [300, -300, 50]
ylist = [0, 0, 200]
names = ['1', '2', '3']
qlist = [1, -20, 10]
ww = Wells(names, xlist, ylist, r_e_m=500)
ww.set_qliq_sm3day(qlist)
xp = np.linspace(-500, 500, 500)
yp = np.linspace(-500, 500, 500)
xv, yv = np.meshgrid(xp, yp)
p_mesh = ww.calc_p_ss_atma(xv, yv)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))
# построим отображение в виде контурной карты
ax1.contour(xp, yp, p_mesh, levels=50)
for x,y,n in zip(xlist, ylist, names):
ax1.text(x+10, y+10, n)
ax1.scatter(xlist, ylist, s =100)
ax2.streamplot(xp, yp, -np.gradient(p_mesh,axis=1),
-np.gradient(p_mesh,axis=0),
density=2)
for x,y,n in zip(xlist, ylist, names):
ax2.text(x+10, y+10, n)
ax2.scatter(xlist, ylist, s =100)
plt.show()
3.9.4 Расчет поля давлений для нескольких скважин на стационарном режиме для заданных забойных давлений
Показать код
xlist = [300, -300, 50]
ylist = [0, 0, 200]
names = ['1', '2', '3']
qlist = [10, -20, 10]
plist = [50, 150, 50]
ww = Wells(names, xlist, ylist, r_e_m=500)
ww.estimate_qliq_from_pwf(names, plist)
xp = np.linspace(-500, 500, 500)
yp = np.linspace(-500, 500, 500)
xv, yv = np.meshgrid(xp, yp)
p_mesh = ww.calc_p_ss_atma(xv, yv)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))
# построим отображение в виде контурной карты
ax1.contour(xp, yp, p_mesh, levels=50)
for x,y,n in zip(xlist, ylist, names):
ax1.text(x+10, y+10, n)
ax1.scatter(xlist, ylist, s =100)
ax2.streamplot(xp, yp, -np.gradient(p_mesh,axis=1),
-np.gradient(p_mesh,axis=0),
density=2)
for x,y,n in zip(xlist, ylist, names):
ax2.text(x+10, y+10, n)
ax2.scatter(xlist, ylist, s =100)
plt.show()
3.10 Производительность скважины
Уравнение производительности скважины можно записать в виде
\[ Q = T \Delta P J_D \tag{3.21}\]
где
\(Q\) - дебит жидкости скважины на поверхности, приведенный к стандартным условиям, м\(^3\)/сут. \[ Q = qB \]
\(T\) - параметр зависящий от гидропроводности пласта \[ T=\dfrac{18.42\mu B q }{\ k h} \tag{3.22}\]
\(\Delta P\) - депрессия на пласт, атм \[ \Delta P = \left(P_e-P_w\right) \tag{3.23}\]
\(J_D\) - безразмерный коэффициент продуктивности скважины, \[ J_D = \dfrac{1}{ \left(\ln{\dfrac{r_e}{r_w}} + S\right)} \tag{3.24}\]
Уравнение 3.21 можно интерпретировать следующим образом. Параметр \(T\) отвечает за свойства пласта и флюида на которые трудно повлиять в ходе эксплуатации. Это то, что дала природа в точке где находится скважина. Депрессия \(\Delta P\) – параметр которым можно управлять в ходе эксплуатации регулируя забойное давление. Например за счет установки насоса и задания параметров его работы. На этом параметре должно быть сосредоточено основное внимание при анализе работы скважины. Параметр \(J_D\) – определяет качество соединения скважины с пластом или качество заканчивания. Его мы можем выбирать при строительстве скважины и можем менять в ходе эксплуатации проводя ГТМ, хотя и достаточно большой ценой. Поскольку мы можем влиять на \(J_D\) важно понимать, какое оптимальное значение продуктивности можно достичь на конкретной скважине и как его можно изменить.
Задачей гидродинамических исследований является установление величин \(T\) и \(J_D\), хотя традиционно речь ведется об определении проницаемости \(k\) и скин-фактора \(S\).
3.11 Задания для самостоятельной работы
Для совершенствования навыков работы с python выполните следующие задания:
Постройте график распределения давления в пласте для композитного пласта. В композитном пласте на расстоянии \(r<r_1\) проницаемость равна \(k=k_1\), а для \(r>=r_1\), \(k=k_2\).
Постройте двумерную тепловую или контурную карту распределения давления в пласте для моделей однородного и композитного пласта.
Рассчитайте среднюю величину давления в круговой области дренирования для однородного пласта. Насколько среднее давление в круговой области дренирования будет отличаться от давления на контуре. Чему будет равен коэффициент \(S\) в выражении \(Q=\dfrac{kh}{18.42\mu B} \dfrac{P_{res}-P_{wf}}{ln(\dfrac{r_e}{r_w})+S}\) при использовании вместо давления на контуре среднего давления? Постройте график, на котором будет отображаться распределение давления в зоне дренирования и величина среднего давления (в виде линии).
Для примера с несколькими скважинами имитирующими трещину ГРП рассчитайте дебиты скважин таким образом, чтобы забойное давление на всех скважинах было одинаковым. Постройте графики распределения давления в пласте. Постройте график дебитов вдоль “скважины”.
4 Нестационарные решения
Для установившегося режима фильтрации давление в пласте не меняется. Для псевдо-установившегося режима постоянным остается перепад давления между пластом и забоем. После запуска, остановки или изменения режима работы скважины эти условия не выполняются. Давление в различных точках пласта может меняться по разному. Такой режим называют неустановившимся, а решения его описывающие нестационарными (зависят от времени).
Неустановившиеся решения уравнения фильтрации (transient solutions) представляют значительный практический интерес во многих задачах, включая задачи интерпретации ГДИС. В тоже время они относительно сложны и требуют применения компьютерных алгоритмов. В данном пособие проведение расчетов иллюстрируется с использованием python и макросов для Excel – Unifloc VBA.
4.1 Решение линейного стока (с использованием Ei)
Уравнение фильтрации для радиального потока в линеаризованном виде можно записать в виде
\[ \frac{\partial p}{\partial t} = 0.00036 \dfrac{k}{\phi \mu c_t} \dfrac{1}{r} \left[ \dfrac{\partial}{\partial r} \left( r \dfrac{ \partial p} {\partial r} \right) \right] \tag{4.1}\]
Напомним, здесь
- \(p\) - давление, атм
- \(t\) - время, час
- \(k\) - проницаемость в направлении движения потока, мД
- \(\mu\) - динамическая вязкость, сП
- \(\phi\) - пористость, д.е.
- \(c_t\) - сжимаемость, 1/атм
- \(r\) - расстояние от центра, м
Часто для анализа уравнений неустановившейся фильтрации используются безразмерные переменные. Мы будем использовать переменные в виде:
\[ r_D = \frac{r}{r_w} \tag{4.2}\] \[ t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2} \tag{4.3}\] \[ p_D = \frac{kh}{ 18.42 q_s B \mu} \left( p_i - p \right) \tag{4.4}\]
Здесь использование практические метрически единицы измерения.
- \(r_w\) - радиус скважины, м
- \(r\) - расстояние от центра скважины до точки в пласте, м
- \(q_s\) - дебит скважины на поверхности, приведенный к нормальным условиям м3/сут
- \(\phi\) - пористость, доли единиц
- \(\mu\) - вязкость нефти в пласте, сП
- \(B\) - объемный коэффициент нефти, м3/м3
- \(p_i\) - начальное давление в пласте, атм
- \(p\) - давление на расстоянии \(r\), атм
- \(c_t\) - общая сжимаемость системы в пласте, 1/атм
Использование безразмерных переменных позволяет упростить уравнение фильтрации, которое примет вид
\[ \frac{\partial p_D}{ \partial t_D} = \dfrac{1}{r_D} \left[ \dfrac{\partial}{\partial r_D} \left( r_D \dfrac{ \partial p_D} {\partial r_D} \right) \right] \tag{4.5}\]
Решение уравнения (Уравнение 4.5) - функция безразмерного давления от безразмерных времени и расстояния \(p_D(r_D, t_D)\).
Для решения уравнения фильтрации - линейного дифференциального уравнения в частных производных второго порядка необходимо задать начальные и граничные условия. Самое простое решение можно получить для случая вертикальной скважины бесконечно малого радиуса запускающейся с постоянным дебитом. Условия соответствующие этому случаю можно выразить следующим образом
начальное условие. До запуска скважины в момент времени \(t_D = 0\) давление в пласте равно начальному во всех точках \(p=p_i\) \[ t_D < 0, p_D = 0 \tag{4.6}\]
условие постоянства дебита на скважине - граничное условие на скважине \[ \lim_{r_D \to 0} {r_D \frac{\partial p_D}{\partial r_D}} = -1 \tag{4.7}\]
условие на бесконечном расстоянии возмущения от скважине нет \[ r_D = \infty, p_D = 0 \tag{4.8}\]
В этом случае решение может быть выражено через функцию интегральной экспоненты \[ p_D(r_D,t_D) = - \frac{1}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right) \tag{4.9}\]
где -Ei(-x) - интегральная показательная функция.
Показать код
# импортируем библиотеки для расчетов
# numpy используем для работы с массивами и подготовки данных для построения графиков.
# Также в некоторых функциях используем возможности векторных расчетов numpy
import numpy as np
# matplotlib используем для построения графиков
import matplotlib.pyplot as plt
# scipy.special используем как альтернативный вариант расчета специальных функций
import scipy.special as sc
# для скорости и удобства используем sc.expi
# построим ветки функции для положительных и отрицательных аргументов раздельно
x = np.arange(1e-5,3,0.01)
x1 = np.arange(-3,-1e-5,0.01)
plt.plot(x, sc.expi(x))
plt.plot(x1, sc.expi(x1))
plt.title("График функции Ei")
plt.xlabel("x")
plt.ylabel("Ei(x)")
plt.show()
где \(-Ei(-x)\) - интегральная показательная функция, рисунок 4.1.
\[ Ei(x)=-\int\limits_{x}^{\infty}\frac{e^{-t}}{t}\,\mathrm dt \]
Часто для проведения расчетов, особенно с использованием компьютерных библиотеке расчетов, бывает удобнее пользоваться модифицированной интегральной показательной функцией \(Ei_1(x)\) или \(E_1(x)\) или \(Ei_n(x)\) при \(n=1\). \[ Ei_n(x) = \int\limits_{1}^{\infty}\frac{e^{-tx}}{t^n}\,\mathrm dt \]
График интегральной показательной функции \(Ei_1(x)\) приведен на Рисунок 4.1. Для вещественных положительных \(x\in\mathbb R, x>0\) верно \(E_1(x) = - Ei( -x)\)
Функцию интегральной экспоненты можно представить в виде ряда.
\[ Ei(x)=-\int\limits_{x}^{\infty}\frac{e^{-t}}{t}\,\mathrm dt=\gamma+\operatorname{ln}|-x|+\sum\limits_{n\ge1}\frac{{-x}^n}{n!\cdot n}, \; x\in\mathbb R,\; \]
Из приведенного выражения можно сделать выводы, что для маленьких значений аргумента функция интегральной экспоненты \(E_1(x)\) может быть аппроксимирована логарифмической зависимостью.
Показать код
# зададим логарифмическое распределение точек вблизи нуля для построения графика
x = np.logspace(-10, 1, 100)
plt.rcParams["figure.figsize"] = (5,9)
fig, (ax1, ax2) = plt.subplots(2, 1)
ax1.plot(x, sc.expi(-x))
ax1.plot(x, np.log(x) + 0.57721566481)
#ax1.set_title('Сравнение Ei и ln для обычных координат')
ax2.plot(x, sc.expi(-x))
ax2.plot(x, np.log(x) + 0.57721566481)
ax2.set_xscale('log')
#ax2.set_title('Сравнение Ei и ln для полулогарифмических координат')
plt.show()
\[ E_1(x) = -ln(x) - \gamma \tag{4.10}\]
График сравнения функций \(E_1(x)\) и \(ln(x)\) показан на рисунке 4.2. Видно, что хорошей аппроксимация будет только для маленьких значений аргумента \(x < 0.01\). Но для решения уравнения фильтрации именно эта зона представляет наибольший интерес.
Представление интегральной экспоненты в виде логарифмической аппроксимации удобно на практике, так как логарифм легче вычислять. В большинстве языков программирования и инструментов для проведения расчетов расчет логарифма реализован по умолчанию. А для расчета интегральной экспоненты, часто приходится предпринимать дополнительные шаги, загружать дополнительные библиотки.
Решение уравнения фильтрации для линейного стока с учетом логарифмической аппроксимации можно представить в виде
\[ p_D(r_D,t_D) = \frac{1}{2} \left( ln \left( \dfrac{ t_D }{r_D^2} \right) +0.809 \right) \tag{4.11}\]
при использовании данного уравнения, следует помнить, что приближенное решение применимо при \(\dfrac{r_D^2}{4t_D} < 0.01\)
Решение линейного стока в размерных переменных
\[ p\left(r,t\right)=p_i-\frac{18.42q_sB\mu}{kh}\left(-\frac{1}{2}Ei\left(-\frac{\varphi\mu c_tr^2}{0.00144kt}\right)\right) \tag{4.12}\]
Решение с учетом логарифмической аппроксимации в размерных переменных
\[ p\left(r,t\right)=p_i-\frac{9.21q_sB\mu}{kh}\left(ln{\frac{kt}{\varphi\mu c_tr^2}}-7.12\right) \tag{4.13}\]
верно при \[ \frac{kt}{\varphi\mu c_tr^2}>70000 \tag{4.14}\]
Решения приведены для практических метрических единиц измерения, что можно увидеть по размерному коэффициенту.
Нестационарное решение с учетом скин-фактора будет иметь вид
\[ P(r, t) = P_{i} - \frac {9.21 {q_s} B\mu }{k h}(\ ln\frac {k t}{ \varphi \mu {c_t} {r^2}} -7.12 + 2S) \tag{4.15}\]
4.2 Построение графиков решения в безразмерных координатах
import scipy.special as sc
"""
Решение линейного стока уравнения фильтрации
"""
def pd_ei(td, rd):
"""
Решение линейного стока уравнения фильтрации
rd - безразмерное расстояние
td - безразмерное время
"""
return -1/2*sc.expi(-rd**2 / 4 / td)Показать код
"""
построим графики решения линейного стока
в безразмерных переменных
"""
rd_arr = np.logspace(1, 3, 100)
td_arr = np.logspace(1, 3, 100)
# при построении используем векторный расчет
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))
for td in [1000, 10000, 100000]:
ax1.plot(rd_arr, pd_ei(td, rd_arr), label=f"td = {td}" )
ax1.set_title("Решение для td = {}".format(td))
ax1.set_xlabel("rd")
ax1.set_ylabel("pd")
ax1.legend()
for rd in (1, 5 , 50):
ax2.plot(td_arr, pd_ei(td_arr, rd), label=f"rd = {rd}" )
ax2.set_title("Решение для rd = {}".format(rd))
ax2.set_xlabel("td")
ax2.set_ylabel("pd")
ax2.legend()
plt.show()
4.3 Построение графиков решения в размерных координатах
Для расчетов в размерных переменных необходимо определить функции перевода размерных координаты в безразмерные и обратно. В целом расчеты станут более громоздкие.
"""
определим функции для перевода размерных переменных в безразмерные и обратно
пригодится для построения графиков и ведения расчетов
"""
def r_from_rd_m(rd, rw_m=0.1):
"""
перевод безразмерного расстояния в размерное
rd - безразмерное расстояние
rw_m - радиус скважины, м
"""
return rd*rw_m
def rd_from_r(r_m, rw_m=0.1):
"""
перевод размерного расстояния в безразмерное
r_m - размерное расстояние, м
rw_m - радиус скважины, м
"""
return r_m/rw_m
def t_from_td_hr(td, k_mD=10, phi=0.2, mu_cP=1, ct_1atm=1e-5, rw_m=0.1):
"""
перевод безразмерного времени в размерное, результат в часах
td - безразмерное время
k_mD - проницаемость пласта, мД
phi - пористость, доли единиц
mu_cP - динамическая вязкость флюида, сП
ct_1atm - общая сжимаемость, 1/атм
rw_m - радиус скважины, м
"""
return td * phi * mu_cP * ct_1atm * rw_m * rw_m / k_mD / 0.00036
def td_from_t(t_hr, k_mD=10, phi=0.2, mu_cP=1, ct_1atm=1e-5, rw_m=0.1):
"""
перевод размерного времени в безразмерное
t_hr - размерное время, час
k_mD - проницаемость пласта, мД
phi - пористость, доли единиц
mu_cP - динамическая вязкость флюида, сП
ct_1atm - общая сжимаемость, 1/атм
rw_m - радиус скважины, м
"""
return 0.00036 * t_hr * k_mD / (phi * mu_cP * ct_1atm * rw_m * rw_m)
def p_from_pd_atma(pd,
k_mD=10, h_m=10,
q_sm3day=20, b_m3m3=1.2, mu_cP=1, pi_atma=250):
"""
перевод безразмерного давления в размерное, результат
в абсолютных атмосферах
pd - безразмерное давление
k_mD - проницаемость пласта, мД
h_m - мощность пласта, м
q_sm3day - дебит на поверхности, м3/сут в с.у.
fvf_m3m3 - объемный коэффициент нефти, м3/м3
mu_cP - динамическая вязкость флюида, сП
pi_atma - начальное давление, абсолютные атм.
"""
return pi_atma - pd * 18.42 * q_sm3day * b_m3m3 * mu_cP / k_mD / h_m
def pd_from_p(p_atma,
k_mD=10, h_m=10,
q_sm3day=20, b_m3m3=1.2, mu_cP=1, pi_atma=250):
"""
перевод размерного давления в безразмерное
p_atma - давление
k_mD - проницаемость пласта, мД
h_m - мощность пласта, м
q_sm3day - дебит на поверхности, м3/сут в с.у.
fvf_m3m3 - объемный коэффициент нефти, м3/м3
mu_cP - динамическая вязкость флюида, сП
pi_atma - начальное давление, абсолютные атм.
"""
return (pi_atma - p_atma) / (18.42*q_sm3day*b_m3m3*mu_cP) * k_mD * h_m Показать код
"""
Построим графики распределения давления и изменения давления
в размерных координатах
"""
# исходные данные для построения графиков
h_m=10
q_sm3day=20
b_m3m3=1.2
mu_cP=1
pi_atma=250
rw_m=0.1
k_mD=10,
phi=0.2
ct_1atm=1e-5
r_arr = np.logspace(0.1, 3, 100)
rd_arr = rd_from_r(r_arr, rw_m=rw_m)
t_arr = np.linspace(1, 100, 100)
td_arr = td_from_t(t_arr,
k_mD=k_mD,
phi=phi,
mu_cP=mu_cP,
ct_1atm=ct_1atm,
rw_m=rw_m)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))
for td in td_from_t(np.array([1, 10, 50]),
k_mD=k_mD,
phi=phi,
mu_cP=mu_cP,
ct_1atm=ct_1atm,
rw_m=rw_m):
ax1.plot(r_arr,
p_from_pd_atma(pd=pd_ei(td,rd_arr),
h_m=h_m,
q_sm3day=q_sm3day,
b_m3m3=b_m3m3,
mu_cP=mu_cP,
pi_atma=pi_atma),
label=f"t = {t_from_td_hr(td, k_mD=k_mD, phi=phi, mu_cP=mu_cP, ct_1atm=ct_1atm, rw_m=rw_m)[0]:.2f} час")
ax1.set_title("Воронка депрессии")
ax1.set_xlabel("r, m")
ax1.set_ylabel("p, atm")
ax1.legend()
for rd in rd_from_r(np.array([rw_m, 10, 50])):
ax2.plot(t_arr,
p_from_pd_atma(pd=pd_ei(td_arr, rd),
h_m=h_m,
q_sm3day=q_sm3day,
b_m3m3=b_m3m3,
mu_cP=mu_cP,
pi_atma=pi_atma),
label=f"r = {r_from_rd_m(rd, rw_m=rw_m)}" )
ax2.set_title("Изменение давления на расстоянии")
ax2.set_xlabel("t, час")
ax2.set_ylabel("p, atm")
ax2.legend()
plt.show()
4.4 Использование библиотеки sympy для построение нестационарного решения
Иногда для манипуляции со сложными выражениями удобно бывает использовать sympy
import sympy as sp
# объявим символьные переменные
r_d, p_d, t_d, gamma = sp.symbols('r_d p_d t_d gamma')
r, p, t, q, b, mu, phi, ct, k, h, p_i, s, r_w = sp.symbols('r p t q b mu phi ct k h p_i s r_w')
# зададим исходное уравнение
eq1 = sp.Eq(p_d, -1/2*sp.Ei(-r_d**2/4 / t_d))
eq1\(\displaystyle p_{d} = - 0.5 \operatorname{Ei}{\left(- \frac{r_{d}^{2}}{4 t_{d}} \right)}\)
решение с интегральной экспонентой можно задать используя scipy.special или сгенерировать из решения sympy
pd_ei_ = sp.lambdify([t_d, r_d],eq1.rhs, modules=["scipy","numpy"])
# сравним решения заданные явно и сгенерированные из решения `sympy`
print(pd_ei(1000,1))
print(pd_ei_(1000,1))
print(pd_ei(1000,1)==pd_ei_(1000,1))3.8585419797881815
3.8585419797881815
TrueПоказать код
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))
rd_arr = np.logspace(0,2,100)
for td in [1, 10, 100, 1000, 10000]:
ax1.plot(rd_arr, pd_ei_(td, rd_arr))
ax1.plot(rd_arr, pd_ei(td, rd_arr))
rd_arr = np.logspace(0,2,100)
for td in [1, 10, 100, 1000, 10000]:
ax2.plot(rd_arr, pd_ei_(td, rd_arr))
ax2.plot(rd_arr, pd_ei(td, rd_arr))
ax2.set_xscale('log')
Связь размерных и безразмерных переменных задается в виде
\[ r_D = \frac{r}{r_w} \] \[ t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2} \] \[ p_D = \frac{kh}{ 18.42 q_s B \mu} \left( p_i - p \right) \]
Здесь использование единицы измерения СИ.
- \(r_w\) - радиус скважины, м
- \(r\) - расстояние от центра скважины до точки в пласте, м
- \(q_s\) - дебит скважины на поверхности, приведенный к нормальным условиям м3/сут
- \(\phi\) - пористость, доли единиц
- \(\mu\) - вязкость нефти в пласте, сП
- \(B\) - объемный коэффициент нефти, м3/м3
- \(p_i\) - начальное давление в пласте, атм
- \(p\) - давление на расстоянии \(r\), атм
- \(c_t\) - общая сжимаемость системы в пласте, 1/атм
# подставим вместо безразмерных переменных размерные
eq2 = eq1.subs(r_d, r/r_w)
eq2 = eq2.subs(t_d, 0.00036*k*t/(phi*mu*ct*r_w**2))
eq2 = eq2.subs(p_d, k*h/(18.42*q*b*mu)*(p_i-p))
# решим полученное уравнение относительно перепада давления
eq3 = sp.solve(eq2,p_i-p)
# выведем первое решение
eq3[0]\(\displaystyle - \frac{9.21 b \mu q \operatorname{Ei}{\left(- \frac{694.444444444444 ct \mu \phi r^{2}}{k t} \right)}}{h k}\)
Решение в размерных переменных можно записать как \[ p\left(r,t\right)=p_i-\frac{18.42q_sB\mu}{kh}\left(-\frac{1}{2} Ei \left(-\frac{\varphi\mu c_tr^2}{0.00144kt}\right)\right) \]
# сгенерируем решение
pd_ei_full_ = sp.lambdify([t, r, k, h, mu, ct, phi, q, b, p_i],p_i-eq3[0],
modules=["scipy","numpy"])
# нарисуем воронку депрессии используя сгенерированное решение
r_arr = np.logspace(-1,2,100)
plt.rcParams["figure.figsize"] = (8,5)
for tt in [0.1, 1, 10]:
plt.plot(r_arr,
pd_ei_full_(tt, r_arr,
k=10, h=10, mu=1, ct=1e-4, phi=0.2, q=10, b=1, p_i=250))
plt.plot(-r_arr,
pd_ei_full_(tt, r_arr,
k=10, h=10, mu=1, ct=1e-4, phi=0.2, q=10, b=1, p_i=250))
Решение с интегральной экспонентой может быть заменено приблеженным решением с использованием логарифма
\[ p_D(r_D,t_D) = - \frac{1}{2} \ln \left( \dfrac{ r_D^2}{4t_d} \right) - \frac{1}{2}\gamma \]
где \(\gamma = 0.57721566481\) - константа Эйлера
на графике от времени в полулогарифмических координатах логарифмическое приближение выглядит как кривая с наклоном \(0.5\)
упростим последнее выражение с использованием sympy
# объявим символьные переменные
r_d, p_d, t_d, gamma = sp.symbols('r_d p_d t_d gamma')
r, p, t, q, b, mu, phi = sp.symbols('r p t q b mu phi')
ct, k, h, p_i, s, r_w = sp.symbols('ct k h p_i s r_w')
# зададим исходное уравнение
eq = sp.Eq(p_d, -1/2*sp.ln(r_d**2/4 / t_d) -1/2*gamma)
print('Исходное уравнение - логарифмическое приближение решения линейного стока')
display(eq)Исходное уравнение - логарифмическое приближение решения линейного стока\(\displaystyle p_{d} = - 0.5 \gamma - 0.5 \log{\left(\frac{r_{d}^{2}}{4 t_{d}} \right)}\)
# подставим значение для gamma
eq = sp.simplify(sp.expand(eq.subs(gamma, 0.57721566)))
eq\(\displaystyle p_{d} = 0.404539350559945 - 0.5 \log{\left(\frac{r_{d}^{2}}{t_{d}} \right)}\)
pd_ei_ln_ = sp.lambdify([t_d, r_d],eq.rhs, modules=["scipy","numpy"]) Показать код
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))
rd_arr = np.logspace(0,2,100)
for td_ in [1, 10, 100, 1000, 10000]:
ax1.plot(rd_arr, pd_ei_(td_, rd_arr))
ax1.plot(rd_arr, pd_ei_ln_(td_, rd_arr))
rd_arr = np.logspace(0,2,100)
for td_ in [1, 10, 100, 1000, 10000]:
ax2.plot(rd_arr, pd_ei_(td_, rd_arr))
ax2.plot(rd_arr, pd_ei_ln_(td_, rd_arr))
ax2.set_xscale('log')
# подставим вместо безразмерных переменных размерные
eq = eq.subs(r_d, r/r_w)
eq = eq.subs(t_d, 0.00036*k*t/(phi*mu*ct*r_w**2))
eq = eq.subs(p_d, k*h/(18.4207*q*b*mu)*(p_i-p))
# решим полученное уравнение относительно перепада давления
eq1 = sp.solve(eq,p_i-p)
# выведем первое решение
eq1[0]\(\displaystyle \frac{9.21035 \cdot 10^{-14} b \mu q \left(80907870111989.0 - 100000000000000.0 \log{\left(\frac{2777.77777777778 ct \mu \phi r^{2}}{k t} \right)}\right)}{h k}\)
упростим решение
eq5 = sp.simplify(sp.expand(eq1[0]))
display(eq5)\(\displaystyle \frac{b \mu q \left(- 9.21035 \log{\left(\frac{ct \mu \phi r^{2}}{k t} \right)} - 65.5807113866197\right)}{h k}\)
4.5 Задания для самостоятельной работы
Для совершенствования навыков работы с python выполните следующие задания:
- Постройте графики зависимости забойного давления от времени в полулогарифмических координатах
- Сравните графики распределения давления вокруг скважины с использованием стационарного решения и решения линейного стока
- Постройте график распределения давления в пласте для композитного пласта. В композитном пласте на расстоянии \(r<r_1\) проницаемость равна \(k=k_1\), а для \(r>=r_1\), \(k=k_2\). Как будет меняться воронка депрессии в таком пласте со временем?
- Постройте двумерную тепловую карту распределения давления в пласте для моделей однородного пласта и композитного пласта.
- Постройте график зависимости радиуса влияния скважины от времени.
- Постройте решение для линейного источника на плоскости проинтегрировав точечное решение по координате. Постройте поле давления для такого источника. Считайте дебит линейного источника известным.
5 Расчет кривой восстановления давления
Один из самых простых примеров применения суперпозиции. Предполагаем, что добывающая скважина в однородном изотропном пласте запускается в момент времени t=0 и работает t_p_hr часов, после чего останавливается. После остановки скважины забойное давление растет - и мы получим кривую восстановления давления.
Пусть решение задачи запуска скважины (падения давления) будет \(P_D(t_D, r_D)\). Тогда решение для изменения давления при запуске и последующей остановки скважины можно представить в виде
\[ P_{bu.D}(t_D, t_{prod.D}, r_D) = P_D(t_D) - P_D(t_D-t_{prod.D}, r_D) \cdot \mathcal{H}(t_D-t_{prod.D}) \tag{5.1}\]
где
- \(t_D\) - безразмерное время после запуска скважины,
- \(t_{prod.D}\) - безразмерное время работы скважины после запуска
- \(\mathcal{H}\) - ступенчатая функция Хевисайда (в некоторых книгах обозначается как \(\theta\))
- \(P_D(t_D, r_D)\) - безразмерное давление - решение задачи запуска скважины (падения давления)
- \(P_{bu.D}(t_D, t_{prod.D}, r_D)\) - безразмерное давление- решение задачи запуска скважины и последующей остановки скважины
Для проведения векторных расчетов в python удобно выражение с использованием функции Хевисайда
\[ \mathcal{H} = \begin{cases}0 & x < 0\\1 & x = 0\\1 & x > 0\end{cases} \tag{5.2}\]
Применение функции Хевисайда позволяет избежать в расчетных функциях применение условных операторов в явном виде для отдельных элементов входных массивов. Это потенциально ускоряет расчет.
def pd_build_up(td, td_p, rd):
"""
расчет давления для запуска и последующей остановки скважины
td - время после запуска
td_p - время безразмерное - которое скважина работала до остановки
rd - расстояния от скважины
"""
# применение функции Хевисайда здесь делает расчет корректным
# для входных векторов td
return pd_ei(td, rd) - np.heaviside(td-td_p,1) * pd_ei(td-td_p, rd)Построим график функции изменения забойного давления от времени.
Показать код
t_arr = np.logspace(-10, 2, 1000)
t_prod_hr = 24
k = 10 # проницаемость
q = 30 # дебит
# переведем размерный массив времени в безразмерные величины
# некоторые параметры можно было бы пропустить и оставить значения по умолчанию
# но для полноты приведем их в явном виде
td_arr = td_from_t(t_arr,
k_mD=k, phi=0.2, mu_cP=1, ct_1atm=1e-05, rw_m=0.1)
td_prod = td_from_t(t_prod_hr,
k_mD=k, phi=0.2, mu_cP=1, ct_1atm=1e-05, rw_m=0.1)
print('время работы скважины {:.2f} часа, что соответсвует безразмерному времени {:.2f}'.format(t_prod_hr, td_prod))
# для заданного массива безразмерных времен рассчитаем безразмерные давления
pd_arr = pd_build_up(td_arr, td_prod, rd=1)
# построение графика
plt.rcParams["figure.figsize"] = (8,3)
plt.plot(td_arr, pd_arr)
plt.xlabel('td')
plt.ylabel('pd')
plt.title('Решение в безразмерных координатах')
plt.show()время работы скважины 24.00 часа, что соответсвует безразмерному времени 4320000.00
Показать код
# переведем безразмерные координаты в размерные
p_arr = p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1, pi_atma=250)
plt.rcParams["figure.figsize"] = (8,3)
plt.plot(t_arr, p_arr)
plt.xlabel('t, hr')
plt.ylabel('p, atma')
plt.title('Решение в размерных координатах')
plt.show()
5.1 Построение графиков распределения давления в пласте при восстановлении давления
Показать код
r_arr = np.logspace(1, 3, 100)
t_arr = np.logspace(-1, 2, 10)
t_arr[0] = 0.001
t_prod = 100
tv, rv = np.meshgrid(td_from_t(t_arr),rd_from_r(r_arr))
pd_arr =pd_ei(tv, rd=rv)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[9,4])
fig.suptitle('Изменение давления в пласте при запуске скважины')
ax1.plot(r_arr,
p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax1.set_label(t_arr)
ax1.set_xlabel('r, m')
ax1.set_ylabel('$p_{wf}, atma$')
ax2.plot(r_arr,
p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax2.set_xscale('log')
ax2.set_xlabel('r, m')
ax2.set_ylabel('$p_{wf}, atma$')
plt.show()
Показать код
pd_arr_bu =pd_ei(tv + td_from_t(t_prod), rd=rv) - pd_ei(tv, rd=rv)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[9,4])
fig.suptitle('Изменение давления в пласте при остановке скважины')
ax1.plot(r_arr,
p_from_pd_atma(pd_arr_bu,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax1.set_label(t_arr)
ax1.set_xlabel('r, m')
ax1.set_ylabel('$p_{wf}, atma$')
ax2.plot(r_arr,
p_from_pd_atma(pd_arr_bu,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax2.set_xscale('log')
ax2.set_xlabel('r, m')
ax2.set_ylabel('$p_{wf}, atma$')
plt.show()
Показать код
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[9,4])
fig.suptitle('Изменение давления в пласте при запуске и остановке скважины')
ax1.plot(r_arr,
p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax1.plot(r_arr,
p_from_pd_atma(pd_arr_bu,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax1.set_label(t_arr)
ax1.set_xlabel('r, m')
ax1.set_ylabel('$p_{wf}, atma$')
ax2.plot(r_arr,
p_from_pd_atma(pd_arr,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax2.plot(r_arr,
p_from_pd_atma(pd_arr_bu,
k_mD=k, q_sm3day=q, h_m=10, b_m3m3=1.2, mu_cP=1,
pi_atma=250))
ax2.set_xscale('log')
ax2.set_xlabel('r, m')
ax2.set_ylabel('$p_{wf}, atma$')
plt.show()
6 Решение для произвольной истории дебитов
6.1 Cтупенчатое изменение дебита
Для расчета изменения давления при переменном дебите введем произвольное референсное значение дебита \(q_{ref}\) (например первое не нулевое значение дебита при запуске скважины). Используем это значение для определения безразмерного давления. \[ p_D = \frac{kh}{ 18.42 q_{ref} B \mu} \left( p_i - p \right) \]
и безразмерного дебита
\[ q_D = \frac{q}{q_{ref}} \]
Тогда, используя принцип суперпозиции, можем выписать выражение для изменения давления на скважине и вокруг нее для произвольного момента времени
\[ P_{mr.D}(t_D, r_D) = \sum_i \left[ q_{D(i)}-q_{D(i-1)} \right] \cdot p_D\left(t_D-t_{D(i)}, r_D\right)\cdot \mathcal{H}(t_D-t_{D(i)}) \tag{6.1}\]
где
- \(i\) - индекс значения дебита в таблице изменения дебитов
- \(q_{D(i)}\) - безразмерный дебит с номером \(i\), который стартует в момент времени \(t_i\). Для первого момента времени \(i\) дебит следующий перед ним считается равным нулю
- \(t_{D(i)}\) - безразмерный момент времени - включения дебита с номером \(i\)
- \(t_{D}\) - безразмерный момент времени для которого проводится расчет
- \(\mathcal{H}\) - ступенчатая функция Хевисайда (Уравнение 5.2)
- \(p_D\left(t\right)\) - зависимость безразмерного давление от времени - решение задачи запуска скважины с постоянным единичным дебитом
- \(P_{mr.D}\) - безразмерное давление \(P_{mr.D}(t_D, r_D)\) учитывающее историю изменения дебитов скважины
# создадим историю изменения дебитов
t_history = np.array([ 0., 2., 24. ], dtype=np.float64)
q_history = np.array([10., 5., 0.], dtype=np.float64)
# массивы должны быть одной и той же длины
def pd_superposition(td, td_hist, qd_hist, rd):
"""
расчет безразмерного давления для последовательности
безразмерных дебитов
td - время расчета после запуска, безразмерное
td_hist - массив времен изменения режимов работы скважин, безразмерное
qd_hist - массив дебитов установленных после изменения режима работы,
безразмерное
"""
# принудительно добавим нули во входные массивы,
# чтобы учесть запуск скважины
qdh = np.hstack([0, qd_hist])
tdh = np.hstack([0, td_hist])
# построим дебиты виртуальных скважин -
# разности реальных дебитов при переключении
delta_qd = np.hstack([0, np.diff(qdh)])
# векторная магия - время как вектор и переключения дебитов тоже
# организуем сумму по временам, внутри сумма по переключениям
# делаем при помощи расчета meshgrid и поиска накопленных сумм
qd_v, td_v =np.meshgrid(delta_qd, td)
# используем куммулятивную сумму numpy для того что суммировать
dpd = np.cumsum(qd_v * pd_ei((td_v - tdh),
rd) * np.heaviside((td_v - tdh), 1), 1)
# последний столбец - полная сумма, которая нужна в качестве результата
return dpd[:,-1]
def q_superposition(t, t_hist, q_hist):
"""
расчет давления для запуска и последующей остановки скважины
t_hr - время после запуска в часах
t_hist_hr - массив времен изменения режимов работы скважин
q_hist_sm3day - массив дебитов установленных после изменения режима работы
"""
def interpolate_constant(x, xp, yp):
indices = np.searchsorted(xp, x, side='right')
y = np.concatenate(([0], yp))
return y[indices]
q=[]
for ti in t:
q.append(interpolate_constant(ti, t_hist, q_hist))
return qПоказать код
td_arr = np.linspace(1e-3, 9e4, 1000)
td_history = np.array([ 0., 2e4, 5e4 ], dtype=np.float64)
qd_history = np.array([10., 5., 0.], dtype=np.float64)
plt.rcParams["figure.figsize"] = (9,5)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(td_arr,
q_superposition(td_arr, td_history, qd_history))
ax2.plot(td_arr,
pd_superposition(td_arr, td_history, qd_history, rd=1),
color='red')
ax2.set_xlabel('td')
ax1.set_ylabel('qd')
ax2.set_ylabel('pd')
plt.show()
def p_superposition_atma(t_hr, t_hist_hr, q_hist_sm3day,
k_mD=10, h_m=10, b_m3m3=1.2, mu_cP=1, pi_atma=250,
phi=0.2, ct_1atm=1e-05, rw_m=0.1):
"""
расчет давления для запуска и последующей остановки скважины
t_hr - время после запуска в часах
t_hist_hr - массив времен изменения режимов работы скважин
q_hist_sm3day - массив дебитов установленных после изменения режима работы
k_mD=10 - проницаемость, мД,
h_m=10 - мощность пласта, м,
b_m3m3=1.2 - объемный коэффициент, м3/м3,
mu_cP=1 - вязкость нефти, сП,
pi_atma=250 - начальное давление, атм,
phi=0.2 - пористость, доли единиц,
ct_1atm=1e-05 - общая сжимаемость, 1/атм,
rw_m=0.1 - радиус скважины
"""
q_ref=1.
td = td_from_t(t_hr, k_mD=k_mD, phi=phi, mu_cP=mu_cP,
ct_1atm=ct_1atm, rw_m=rw_m)
td_hist = td_from_t(t_hist_hr, k_mD=k_mD, phi=phi, mu_cP=mu_cP,
ct_1atm=ct_1atm, rw_m=rw_m)
return p_from_pd_atma(pd_superposition(td, td_hist, q_hist_sm3day / q_ref, 1),
k_mD=10, h_m=10,
q_sm3day=q_ref, b_m3m3=1.2, mu_cP=1, pi_atma=250)Показать код
t_arr = np.arange(1e-6, 50, 1e-2)
plt.rcParams["figure.figsize"] = (9,5)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(t_arr, q_superposition(t_arr, t_history, q_history))
ax2.plot(t_arr, p_superposition_atma(t_arr, t_history, q_history), color='red')
ax2.set_xlabel('t, hr')
ax1.set_ylabel('q, m3/day')
ax2.set_ylabel('p, atma')
plt.show()
Показать код
t_history = np.array([0, 2, 10, 24, 24.1,24.2,24.3,24.4,24.5 ])
q_history = np.array([10, 11, 12, 11, 10, 9, 8, 7, 6])
plt.rcParams["figure.figsize"] = (9,5)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(t_arr, q_superposition(t_arr, t_history, q_history))
ax2.plot(t_arr, p_superposition_atma(t_arr, t_history, q_history, 1), color='red')
ax2.set_xlabel('t, hr')
ax1.set_ylabel('q, m3/day')
ax2.set_ylabel('p, atma')
plt.show()
6.2 Линейное изменение дебита
Для линейно меняющегося дебита \(dq_D t_D\) решение можно представить в виде интеграла
\[ p_D = \int_0^{t_D}{- \frac{dq_D}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right) dt_D} \tag{6.2}\]
Выражение (6.2) аналитически интегрируется, что можно легко показать используя библиотеку sympy.
rd, td, qd = sp.symbols('r_D t_D q_D')
# запишем решение в символьном виде для постоянного дебита
eq = -qd/2 * sp.Ei(-rd **2 / 4 / td)
display(eq)\(\displaystyle - \frac{q_{D} \operatorname{Ei}{\left(- \frac{r_{D}^{2}}{4 t_{D}} \right)}}{2}\)
# проинтегрируем решение по времени
eq1 = sp.integrate(eq, td)
display(eq1)\(\displaystyle - \frac{q_{D} \left(\frac{r_{D}^{2} \operatorname{Ei}{\left(- \frac{r_{D}^{2}}{4 t_{D}} \right)}}{4} + t_{D} \operatorname{Ei}{\left(- \frac{r_{D}^{2}}{4 t_{D}} \right)} + t_{D} e^{- \frac{r_{D}^{2}}{4 t_{D}}}\right)}{2}\)
Для линейной зависимости дебита от времени \[ Q_D = dQ_D \cdot t_D \tag{6.3}\]
можно получить выражение
\[ p_D(r_D,t_D, dQ_D) =-\frac{dQ_D t_D }{2} \left[ \left( 1+ \frac{r_D^2}{4 t_D} \right) Ei \left(- \dfrac{r_D^2}{4t_D} \right) + e^{-\dfrac{r_D^2}{4t_D}} \right] \tag{6.4}\]
где \(dQ_D\) - скорость изменения дебита.
Для таблично заданных дебитов и времен можно оценить
\[ dQ_{D(i)} = \dfrac{Q_{D(i)}-Q_{D(i-1)}}{t_{D(i)} - t_{D(i-1)} } \]
Cравните формулу (6.4) с формулой (9.68) в книге Щелкачева “Основы неустановившейся фильтрации” (Щелкачев 1995)
Щелкачев, В. Н. 1995. Основы и приложения теории неустановившейся фильтрации. Часть 1. Москва: Нефть и газ.
Тогда, используя принцип суперпозиции, можем выписать выражение для изменения давления на скважине и вокруг нее для произвольного момента времени
\[ P_{mr.D}(t_D, r_D) = \sum_i p_D\left(t_D-t_{D(i)}, r_D, dQ_{D(i+1)} - dQ_{D(i)}\right)\cdot \mathcal{H}(t_D-t_{D(i)}) \tag{6.5}\]
где
- \(i\) - индекс значения дебита в таблице изменения дебитов
- \(dQ_{D(i)}\) - изменение безразмерного дебита относительно безразмерного времени (14.4)
- \(t_{D(i)}\) - безразмерный момент времени - включения дебита с номером \(i\)
- \(t_{D}\) - безразмерный момент времени для которого проводится расчет
- \(\mathcal{H}\) - ступенчатая функция Хевисайда
- \(p_D\left(t\right)\) - зависимость безразмерного давление от времени - решение задачи запуска скважины с постоянным единичным дебитом
- $P_{mr.D} $ - безразмерное давление \(P_{mr.D}(t_D, r_D)\) учитывающее историю изменения дебитов скважины
следует обратить внимание, при суперпозиции скорость изменения дебита вычисляется как \(dQ_{D(i+1)} - dQ_{D(i)}\). при реализации расчета необходимо предусмотреть, чтобы для первого и последнего шага расчет прошел корректно. Для этого можно, например, добавить к массивам дебитов и времени дополнительный значения в начале и в конце массивов соответствующие постоянным значениям дебита. Также надо учитывать, что в приведенном выражении массивы должны начинаться со значений \(Q_D=0\)
# Решение линейного стока уравнения фильтрации
def pd_lin_ei(td, rd=1, dqd_dtd=1):
"""
Решение линейного стока уравнения фильтрации
rd - безразмерное расстояние
td - безразмерное время
"""
# при расчете убедимся, что td=0 не повлияет на расчет,
# даже если td массив и нулевой только один элемент
td = np.array(td, dtype = float)
pd = (1 + rd**2/4/td) * (-sc.expi(-rd**2/4/td)) - np.exp(-rd**2/4/td)
return dqd_dtd * td * pd / 2
def pd_superposition_lin(td, td_hist, qd_hist):
"""
расчет безразмерного давления для последовательности
безразмерных дебитов
- td - время расчета после запуска, безразмерное
- td_hist - массив времен изменения режимов работы, безразмерное
- qd_hist - массив дебитов после изменения режима работы, безразмерное
"""
# принудительно добавим нули во входные массивы,
# чтобы учесть запуск скважины
qdh = np.hstack([qd_hist])
tdh = np.hstack([td_hist])
# построим дебиты виртуальных скважин
# разности реальных дебитов при переключении
delta_qd = np.hstack([np.diff(qdh),0])
delta_td = np.hstack([np.diff(tdh),1])
dq_dt = delta_qd / delta_td
dq_dt = np.diff(np.hstack([0, delta_qd / delta_td]))
# референсный безразмерный дебит это 1
# векторная магия - время как вектор и переключения дебитов тоже
# организуем сумму по временам, внутри сумма по переключениям
# делаем при помощи расчета meshgrid и поиска накопленных сумм
qd_v, td_v =np.meshgrid(delta_qd, td)
dpd = np.cumsum(pd_lin_ei((td_v - tdh),
dqd_dtd=dq_dt) * np.heaviside((td_v - tdh), 1),1 )
return dpd[:,-1]
def p_superposition_lin_atma(t_hr, t_hist_hr, q_hist_sm3day,
k_mD=10, h_m=10, b_m3m3=1.2, mu_cP=1, pi_atma=250,
phi=0.2, ct_1atm=1e-05, rw_m=0.1):
"""
расчет давления для запуска и последующей остановки скважины
t_hr - время после запуска в часах
t_hist_hr - массив времен изменения режимов работы скважин
q_hist_sm3day - массив дебитов установленных после изменения режима работы
k_mD=10 - проницаемость, мД,
h_m=10 - мощность пласта, м,
b_m3m3=1.2 - объемный коэффициент, м3/м3,
mu_cP=1 - вязкость нефти, сП,
pi_atma=250 - начальное давление, атм,
phi=0.2 - пористость, доли единиц,
ct_1atm=1e-05 - общая сжимаемость, 1/атм,
rw_m=0.1 - радиус скважины
"""
q_ref = 1.
td = td_from_t(t_hr,
k_mD=k_mD, phi=phi, mu_cP=mu_cP, ct_1atm=ct_1atm, rw_m=rw_m)
td_hist = td_from_t(t_hist_hr,
k_mD=k_mD, phi=phi, mu_cP=mu_cP,
ct_1atm=ct_1atm, rw_m=rw_m)
return p_from_pd_atma(pd_superposition_lin(td, td_hist, q_hist_sm3day / q_ref),
k_mD=10, h_m=10,
q_sm3day=q_ref, b_m3m3=1.2, mu_cP=1, pi_atma=250)Показать код
td_arr = np.linspace(1e-3, 2000, 2000)
td_history = np.array([0., 1, 50, 100, 300, 500 ], dtype=np.float64)
qd_history = np.array([0., 3, 1.4, 1.0, 1.01, 1], dtype=np.float64)
plt.rcParams["figure.figsize"] = (9,6)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(td_arr, q_superposition(td_arr, td_history, qd_history))
ax1.plot(td_arr, np.interp(td_arr, td_history, qd_history))
ax2.plot(td_arr,
pd_superposition_lin(td_arr, td_history, qd_history),
color='red')
ax2.plot(td_arr,
pd_superposition(td_arr, td_history, qd_history, rd=1),
color='blue')
ax2.set_xlabel('td')
ax1.set_ylabel('qd')
ax2.set_ylabel('pd')
plt.show()
Показать код
t_history = np.array([0, .02, 10, 24, 24.1,24.2,24.3,24.4,24.5 ])
q_history = np.array([0, 11, 12, 12, 10, 9, 8, 7, 6])
plt.rcParams["figure.figsize"] = (9,6)
fig, (ax1, ax2) = plt.subplots(2,1)
ax1.plot(t_arr, q_superposition(t_arr, t_history, q_history))
ax1.plot(t_arr, np.interp(t_arr, t_history, q_history))
ax2.plot(t_arr,
p_superposition_atma(t_arr, t_history, q_history),
color='red')
ax2.plot(t_arr,
p_superposition_lin_atma(t_arr, t_history, q_history),
color='green')
ax2.set_xlabel('t, hr')
ax1.set_ylabel('q, m3/day')
ax2.set_ylabel('p, atma')
plt.show()
7 Радиус влияния скважины
Нестационарное решение в безразмерных переменных \[ p_D(r_D,t_D) = - \frac{1}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right) \]
где безразмерные переменные введены как
\[ r_D = \frac{r}{r_w} \]
\[ t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2} \]
\[ p_D = \frac{kh}{ 18.42 q_s B \mu} \left( p_i - p \right) \]
Здесь использование единицы измерения СИ.
- \(r_w\) - радиус скважины, м
- \(r\) - расстояние от центра скважины до точки в пласте, м
- \(q_s\) - дебит скважины на поверхности, приведенный к нормальным условиям м3/сут
- \(\phi\) - пористость, доли единиц
- \(\mu\) - вязкость нефти в пласте, сП
- \(B\) - объемный коэффициент нефти, м3/м3
- \(p_i\) - начальное давление в пласте, атм
- \(p\) - давление на расстоянии \(r\), атм
- \(c_t\) - общая сжимаемость системы в пласте, 1/атм
Для этих же безразмерных переменных, считая начальное давление равным давлению на контуре можно записать стационарное решение для давления в круговом пласте
\[ p_D = \ln r_{eD} - \ln r_D \]
сравним это решение с логарифмической аппроксимацией (4.11)
\[ p_D(r_D,t_D) = - \frac{q_D}{2} \left[ \ln \left( \dfrac{ r_D^2}{4t_d} \right) +\gamma \right] \tag{7.1}\]
которое можно преобразовать к виду
\[ p_D(r_D,t_D) = - q_D \ln r_D + \frac{q_D}{2} \left[ \ln(4t_D) -\gamma \right] \tag{7.2}\]
сравнивая со стационарным решением можно найти выражение безразмерного радиуса контура в зависимости от безразмерного времени \[ \ln r_{eD} = \frac{1}{2}(\ln(4t_D)-\gamma) \tag{7.3}\]
\[ r_{eD} = \sqrt { 4t_D e^{-\gamma} } \tag{7.4}\]
# оценим значение величины под корнем
print(4*np.exp(-0.57721566481))2.2458379344731085наконец получим \[ r_{eD} = \sqrt {2.2458 t_D} \tag{7.5}\]
это значение называют радиусом влияния скважины. Используя это значение для определенного момента времени можно получить стационарное распределение давления в системе хорошо приближающее решение линейного стока работающего в бесконечном пласте. Можно считать это расстояние за расстояние на которое распространяется влияние скважины.
достижение радиуса влияния внешних границ будет обуславливать начало перехода от неустановившегося режима фильтрации к режиму обусловленному влиянием границ - стационарному для границы постоянного давления или псевдоустановившемуся для границы непротекания.
# Решение линейного стока уравнения фильтрации
def pd_ei(td, rd):
"""
Решение линейного стока уравнения фильтрации
- rd - безразмерное расстояние
- td - безразмерное время
"""
return -1/2*sc.expi(-rd**2 / 4 / td)
def pd_ss(rd, red):
"""
стационарное решение в безразмерных переменных
- rd - безразмерное расстояние
- red - безразмерное расстояние до границы
"""
return np.log(red/rd)Показать код
# зададим точки расстояний для отрисовки графика
rdl = np.logspace(-1, 3 , 100)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[9,4])
# построим первый график в обычных координатах
for td in [10, 100, 1000]:
ax1.plot(rdl, pd_ei(td, rdl))
red = (td*2.2458)**0.5
ax1.plot(rdl, pd_ss(rdl, red))
# построим второй график в полулогарифмических координатах
for td in [10, 100, 1000]:
ax2.plot(rdl, pd_ei(td, rdl))
red = (td*2.2458)**0.5
ax2.plot(rdl, pd_ss(rdl, red))
plt.xscale('log')
plt.show()
8 Преобразования размерных величин с использованием python и sympy
Преобразования размерных величин удобно выполнять с модулем символьных вычислений python - sympy. Преобразования размерностей ключевых величин полезно знать наизусть, хотя всегда можно найти их в таблицах. Значения многих физические константы зашины в модуле scipy.constants, иногда это оказывается удобным, при этом автоматически будет учитываться достаточно большое количество знаков после запятой в константах. Рассмотрим размерности ряда величин широко применяемых в нефтяном инжиниринге.
# импортируем sympy и константы из scipy
import sympy as sp
import scipy.constants as const
# в модуле scipy.constants есть значения общепринятых констант -например значение pi
print(const.pi)3.1415926535897938.1 Задание размерных величин
8.1.1 Объемный расход \(q\)
В СИ измеряется в [м\(^3\)/сек], в практических метрических единицах измеряется в [м\(^3\)/сут], в американских промысловых единицах измеряется в [bbl/day].
- \(1\) [м\(^3\)/сек] = \(543439\) [bbl/day] = \(86400\) [м\(^3\)/сут]
- \(1\) [м\(^3\)/сут] = \(\dfrac{1}{86400}\) [м\(^3\)/сек] \(= 1.157407 \cdot 10^{-5}\) [м\(^3\)/сек]
- \(1\) [bbl/day] = \(0.15898\) [м\(^3\)/сут]
# выведем некоторые переводные коэффициенты для объемных расходов
print(f'Одни [сут] = {24*60*60} = {const.day} [сек]')
print(f'Один [м3/сут] = {1/const.day} [м3/сек]')
print(f'Один баррель в день [bbl/day] = {const.bbl} [м3/сут]')
print(f'Один баррель в день [bbl/day] = {const.bbl/const.day} [м3/сек]')
print(f'Один [м3/сут] = {1/const.day} [м3/сек]')Одни [сут] = 86400 = 86400.0 [сек]
Один [м3/сут] = 1.1574074074074073e-05 [м3/сек]
Один баррель в день [bbl/day] = 0.15898729492799998 [м3/сут]
Один баррель в день [bbl/day] = 1.8401307283333331e-06 [м3/сек]
Один [м3/сут] = 1.1574074074074073e-05 [м3/сек]8.1.2 Проницаемость \(k\)
В СИ измеряется в [м\(^2\)], в практических метрических единицах измеряется в [мД], в американских промысловых единицах измеряется в [mD].
Определение: в пористой среде с проницаемостью в один Дарси для поддержания течения жидкости с динамической вязкостью 1 сП со скоростью фильтрации 1 см/с необходимо поддерживать перепад давления жидкости приблизительно в одну атмосферу на 1 см вдоль направления течения. При использовании физической атмосферы для расчета перепада давления (физическая атмосфера = 101 325 Па) единица проницаемости равняется приблизительно 0.986923 мкм².
В отечественной литературе при определении дарси в качестве величины атмосферы было принято использовать техническую атмосферу (1 кгс/см² = 98 066,5 Па), так что для величины дарси получалось значение приблизительно 1,02 мкм², причём эпизодические случаи использования западного определения дарси специально отмечались [ru.wikipedia.org/wiki/Дарси]. Согласно ГОСТ 26450.2-85 величины 1 Дарси \(= 0.9869⋅10^{−12}\) м².
- \(1\) [м\(^2\)] = \(1.01325 \cdot 10^{15}\) [мД]
- \(1\) [мД] = \(0.986923 \cdot 10^{-15}\) [м\(^2\)]
print(f'Один [мД] = {1e5/const.atm * 1e-15} [м²]')Один [мД] = 9.86923266716013e-16 [м²]8.1.3 Вязкость \(\mu\)
- \(1\) [Па\(\cdot\)с] = \(1000\) [сП]
- \(1\) [сП] = \(10^{-3}\) [Па\(\cdot\)с]
8.1.4 Давление \(p\)
- \(1\) [Па] = \(0.0001450\) [psi] = \(0.00000987\) [атм]
- \(1\) [атм] = \(14.6959\) [psi] = \(101325\) [Па]
AT = 98066.5 # technical atmosphere in Pa, техническая атмосфера в Па
print(f'Один [psi] в [Па] = {const.psi}')
print(f'Один [bar] в [Па] = {const.bar}')
print(f'Один [atm] в [Па] = {const.atm}')
print(f'Один [at] в [Па] = {AT}')
print(f'Один [atm] в [psi] = {const.atm/const.psi}')Один [psi] в [Па] = 6894.757293168361
Один [bar] в [Па] = 100000.0
Один [atm] в [Па] = 101325.0
Один [at] в [Па] = 98066.5
Один [atm] в [psi] = 14.695948775513458.1.5 Расстояние \(x\)
- \(1\) [м] = \(3.28\) [ft]
8.2 Размерный коэффициент для формулы Дюпюи
Используя рассчитанные выше переводные коэффициенты для различных размерных величин рассчитаем переводной коэффициент в формуле Дюпюи
\[ Q = \frac{ 2 \pi kh}{ \mu B} \frac{ \left( p_i - p \right) } {\ln{\dfrac{r_e}{r_w}} +S } \]
# зададим переменные sympy
Q, k, h, mu, B, pres, pwf, re, rw, S, pi = sp.symbols('Q k h mu B p_res p_wf r_e r_w S pi',
real=True,
positive=True)
# определим уравнение
eq = sp.Eq(Q, 2 * pi * k * h / (mu * B) * (pres - pwf) / (sp.ln(re/rw) + S))
display(eq)\(\displaystyle Q = \frac{2 h k \pi \left(p_{res} - p_{wf}\right)}{B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
eq = eq.subs(Q, 1/const.day * Q) # дебит, [м3/сут] в [м3/сек]
f_k = 1e5/const.atm * 1e-15
eq = eq.subs(k, f_k * k) # проницаемость, [мД] в [м2]
eq = eq.subs(mu, 1e-3 * mu) # вязкость, [сП] в [Па сек]
eq = eq.subs(pres, const.atm * pres) # давление [атм] в [Па]
eq = eq.subs(pwf, const.atm * pwf) # давление [атм] в [Па]
eq = eq.subs(pi, const.pi)
display(eq)\(\displaystyle 1.15740740740741 \cdot 10^{-5} Q = \frac{6.20102176874373 \cdot 10^{-12} h k \left(101325.0 p_{res} - 101325.0 p_{wf}\right)}{B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
Решим полученное уравнение относительно Q и упростим средствами sympy
eq1 = sp.simplify(sp.solve(eq,Q)[0])
display(sp.Eq(Q,eq1))\(\displaystyle Q = \frac{0.0542867210540316 h k \left(p_{res} - p_{wf}\right)}{B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
Выделим полученную константу в явном виде и найдем обратную величину - это и будет необходимый нам переводной коэффициент.
f = 1/eq1.args[0]
print(f)18.4207110060064По умолчанию sympy автоматически организует порядок элементов в своих выражениях. Этот порядок может отличаться от привычного - хотя и суть формул при этом не меняется. Применяя некоторые хитрости можно заставить sympy вывести выражения в приемлимом виде.
a = sp.symbols('a')
eq2 = eq1.subs(eq1.args[0],1/a)
with sp.evaluate(False):
display(eq2.subs(a, f))\(\displaystyle \frac{h k \left(p_{res} - p_{wf}\right)}{18.4207110060064 B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
Но иногда результат проще переписать руками в нужном виде. В итоге уравнение Дюпюи в практических метрических единицах измерения примет вид.
\[ Q = \frac{ kh}{18.42 \mu B} \frac{ \left( p_i - p \right) } {\ln{\dfrac{r_e}{r_w}} +S } \]
где
- \(Q\) - дебит скважины на поверхности, приведенный к нормальным условиям, ст. м\(^3\)/сут
- \(\mu\) - вязкость нефти в пласте, сП
- \(B\) - объемный коэффициент нефти, м\(^3\)/м\(^3\)
- \(P_{res}\) - пластовое давление или давление на контуре с радиусом \(r_e\), атма
- \(P_{wf}\) - давление забойное, атма
- \(k\) - проницаемость, мД
- \(h\) - мощность пласта, м
- \(r_e\) - внешний контур дренирования скважины, м
- \(r_w\) - радиус скважины, м
- \(S\) - скин-фактор скважины, безразмерн.